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Abstract 

On  the  Numerical  Solution  of  One-Dimensional 
Integral  and  Differential  Equations 


Harold  Page  Starr,  Jr. 
Yale  University 
1992 


Many  problems  in  mathematical  physics  can  be  formulated  as  one- dimensional  integral 
equations.  Examples  include  problems  in  electrostatics,  crack  problems  in  elastic  bodies, 
and  two-point  boundary  value  problems  for  ordinary  differential  equations.  Since  most  in¬ 
tegral  equations  arising  in  applications  do  not  have  analytic  solutions,  there  is  considerable 
interest  in  the  numerical  solution  of  these  problems.  Unfortunately,  discretization  of  inte¬ 
gral  equation  ads  to  dense  systems  of  linear  algebraic  equations,  and  the  direct  solution 
of  a  dense  linear  system  of  dimension  N  requires  order  O(N^)  arithmetic  operations.  Al¬ 
ternatively,  the  solution  to  the  linear  system  can  be  obtained  using  an  iterative  method 
such  as  the  conjugate  gradient  algorithm  or  the  conjugate  residual  algorithm.  When  the 
condition  number  of  the  linear  system  is  small,  the  amount  of  work  required  is  reduced  to 
O(IV^),  and  can  be  reduced  further  to  0(N)  when  the  iterative  method  is  combined  with 
an  algorithm  such  as  the  fast  multipole  method.  However,  many  problems  (including  those 
formulated  as  first  kind  integral  equations)  yield  ill-conditioned  linear  systems;  for  these 
problems,  the  cost  of  an  iterative  method  is  prohibitive,  even  when  combined  with  an  algo¬ 
rithm  such  as  the  fast  multipole  method.  Recently,  wavelet-like  bases  have  been  developed 
with  the  property  that  integral  operators  in  these  bases  correspond  to  matrices  which  are 
sparse.  When  the  inverses  of  these  integral  operators  also  correspond  to  sparse  matrices, 
the  Schulz  method  becomes  highly  effective  and  produces  an  0{N  -log^  N)  algorithm  for 
solving  a  one-dimensional  integral  equation.  Unfortunately,  for  first  kind  integral  equations 
and  other  problems  of  interest,  the  integral  operators  do  not  have  sparse  inverses  in  these 
wavelet-like  bases. 

This  thesis  is  based  on  the  observation  that  one-dimensional  integral  operators  can 
be  recursively  decomposed  into  sums  of  products  of  operators  of  low  numerical  rank.  A 
complicated  analytical  apparatus  is  then  constructed  which  allows  for  the  direct  solution  of 
an  integral  equation  in  order  0{N)  operations.  The  algorithms  of  this  thesis  permit  the  use 
of  schemes  with  extremely  high  orders  of  convergence,  and  are  quite  insensitive  to  end-point 
singularities.  The  performance  of  the  methods  is  illustrated  with  numerical  examples. 
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Chapter  1 

Introduction 


Many  problems  in  mathematical  physics  can  be  formulated  as  one-dimensional  integral 
equations.  From  an  abstract  viewpoint,  the  advantage  of  the  integral  equation  formulation  is 
that  many  properties  of  the  solution  are  readily  apparent;  from  a  computational  viewpoint, 
there  exist  extremely  stable,  high  order  numerical  methods  for  the  solution  of  integral 
equations.  In  addition,  linear  systems  which  arise  from  discretization  of  second  kind  integral 
equations  are  generally  well-conditioned.  On  the  other  hand,  linear  systems  arising  from 
first  kind  integral  equations  generally  have  condition  numbers  of  at  least  0{N),  v/here  N  is 
the  number  of  points  in  the  discretization. 

Despite  their  advantages,  integral  equations  are  virtually  never  used  as  a  numerical 
tool,  since  their  discretization  leads  to  dense  systems  of  linear  algebraic  equations,  and  the 
solution  of  a  dense  linear  system  of  dimension  N  reqtiires  order  O(N^)  arithmetic  operations, 
with  N  the  number  of  nodes  in  the  discretization.  This  makes  the  use  of  integral  equations 
extremely  unattractive  as  a  numerical  tool,  despite  their  desirable  analytical  properties. 

In  recent  years,  a  number  of  algorithms  has  been  developed  for  the  fast  application  of 
integral  operators  [4],  [21],  the  best  known  of  which  are  the  particle  simulation  algorithms 
developed  by  L.  Greengard  and  V.  Rokhlin.  Each  algorithm  of  this  clciss  exploits  the  special 
structure  of  a  particular  problem  by  combining  interpolation  of  the  function  which  defines 
the  matrix  elements  with  a  divide-and-conquer  strategy,  leading  to  a  scheme  for  applying 
the  integral  operator  to  an  arbitrary  vector  for  a  cost  proportional  to  N  (or,  sometimes, 
N  •  logiV),  where  N  is  the  number  of  elements  in  the  discretization  of  the  domain  of  the 
operator. 

When  such  a  scheme  is  combined  with  a  conjugate  gradient  type  procedure  for  the 
solution  of  the  integral  equation,  the  resulting  algorithm  requires  (asymptotically)  a  finite 
number  of  iterations  to  converge,  leading  to  an  order  0{N)  estimate  for  the  solution  of  the 
original  integral  equation.  Unfortunately,  the  actual  number  of  iterations  required  is  very 
sensitive  to  the  conditioning  of  the  problem  being  solved.  In  many  cases  of  interest,  this 
number  is  prohibitive. 

In  [5],  orthonormal  bases  are  developed  with  the  property  that  integral  operators  in 
these  bases  correspond  to  matrices  which  are  sparse.  When  their  inverses  also  correspond 
to  sparse  matrices,  the  Schulz  method  becomes  highly  effective  and  produces  an  order 
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0(N  -log^  N)  algorithm  for  the  solution  of  a  second  kind  integral  equation.  While,  formally 
speaking,  the  Schultz  technique  is  an  iterative  one,  in  reality  it  behaves  almost  like  a  direct 
algorithm,  since  the  number  of  iterations  it  requires  is  proportional  to  the  logarithm  of 
the  condition  number  k  of  the  matrix  (as  opposed  to  \/k  and  k  for  the  conjugate  gradient 
and  conjugate  residual  techniques,  respectively).  However,  this  method  can  not  be  used  for 
operators  which  do  not  have  sparse  inverses  (for  example,  first  kind  integral  operators). 

In  [22],  a  direct  method  was  developed  for  the  solution  of  second  kind  integral  equations 
resulting  from  two-point  boundary  value  problems  of  second  order  ordinary  differential 
equations.  For  these  problems,  the  integral  operators  can  be  recursively  decomposed  into 
sums  of  products  of  operators  of  low  rank.  A  somewhat  involved  analytical  apparatus  is 
then  constructed  which  allows  for  the  direct  solution  of  the  integral  equation  in  order  0(iV) 
operations,  with  N  the  number  of  nodes  on  the  interval. 

In  this  thesis,  we  construct  0(N)  algorithms  for  the  direct  solution  of  first  and  second 
kind  integral  equations.  We  first  extend  [22]  to  permit  the  fast,  direct  solution  of  two-point 
boundary  problems  of  systems  of  first  order  ordinary  differential  equations.  We  then  extend 
the  observations  of  [5]  to  construct  sparse  representations  of  integral  operators  with  either 
weakly  singular  kernels  or  a  Cauchy  kernel,  and  extend  the  techniques  of  [22]  and  [29]  to 
apply  the  inverse  of  these  operators  using  order  0(N)  arithmetic  operations. 

The  plan  of  this  thesis  is  as  follows:  Chapter  2  describes  the  algorithms  for  two-point 
boundary  value  problems  for  systems  of  ordinary  differential  equations  (this  chapter  has 
been  published  previously  [29]);  Chapter  3  describes  the  algorithms  for  integral  equations 
with  singularities;  Chapter  4  briefly  outlines  some  generalizations,  and  presents  our  conclu¬ 
sions. 
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Two-Point  Boundary  Value 
Problems 


Second  kind  integral  equations  have  been  a  popular  analytical  tool  in  the  study  of  ordi¬ 
nary  differential  equations  for  nearly  a  century.  When  boundary  value  problems  are  being 
considered,  the  integral  equations  which  arise  are  of  the  Fredholm  type.  From  an  abstract 
viewpoint,  the  advantage  of  this  formulation  is  that  many  properties  of  the  solution  are 
readily  apparent;  from  a  computational  viewpoint,  the  linear  systems  which  arise  from  dis¬ 
cretization  are  generally  well-conditioned.  An  ill-behaved  differential  equation  can  often 
be  reduced  to  a  perfectly  tractable  integral  equation  by  means  of  an  appropriate  choice 
of  the  “background”  Green’s  function  (see  Example  2.3  in  Section  2.4  b'^low).  Standard 
finite  difference  and  finite  element  methods,  on  the  other  hand,  which  discretize  the  original 
differential  equation,  encounter  serious  numerical  difficulties  when  the  solution  possesses 
derivatives  of  large  magnitude  (boundary  layers).  A  second  advantage  is  that  there  exist 
extremely  stable,  high  order  numerical  methods  for  the  solution  of  second  kind  Fredholm 
equations,  while  the  order  of  convergence  of  most  practical  schemes  for  the  solution  of 
ordinary  differential  equations  tends  to  be  limited,  even  if  Richardson  extrapolation  and 
deferred  correction  approaches  are  considered. 

Despite  all  these  advantages,  integral  equations  are  virtually  never  used  as  a  numerical 
tool  for  the  solution  of  systems  of  two-point  boundary  value  problems,  since  their  discretiza¬ 
tion  leads  to  dense  systems  of  linear  algebraic  equations,  and  the  solution  of  a  dense  linear 
system  of  dimension  N  ■  n  requires  order  0{N^  ■  n^)  arithmetic  operations,  with  N  the 
number  of  nodes  in  the  discretization,  and  n  the  number  of  equations  in  the  system.  Finite 
difference  and  finite  element  schemes  lead  to  banded  systems  of  linear  algebraic  equations, 
and  the  solution  of  the  latter  requires  order  0{N-n^)  arithmetic  operations.  This  makes  the 
use  of  integral  equations  extremely  unattractive  as  a  numerical  tool,  despite  their  superior 
analytical  properties.  A  similar  difficulty  is  encountered  when  spectral  methods  are  applied 
to  boundary  value  problems.  They  yield  high  order  accuracy,  but  result  in  dense  systems 
of  linear  algebraic  equations. 

Recently.  [22]  presented  a  fast  numerical  algorithm  for  solving  two-point  boundary  value 
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problems  for  second  order  differential  equations.  By  solving  the  problems  as  second  kind 
integral  equations,  one  obtains  the  superior  properties  of  integral  equations  over  differential 
equations.  By  using  the  technique  of  [22],  integral  equations  arising  from  boundary  value 
pioblems  are  solved  in  order  0{N  •  p^)  arithmetic  operations,  with  N  the  number  of  nodes 
on  the  interval  and  p  the  desired  order  of  convergence. 

In  this  chapter,  we  extend  the  results  of  [22]  by  showing  that  integral  equations  arising 
from  two-point  boundary  value  problems  for  systems  of  ordinary  differential  equations  can 
be  solved  in  0(N  ■  p^  •  n^)  arithmetic  operations,  with  n  the  number  of  equations  in  the 
system.  We  in  addition  present  a  Newton  method  for  solving  boundary  value  problems  for 
nonlinear  first  order  systems  in  which  each  Newton  iterate  is  the  solution  of  a  second  kind 
integral  equation. 

The  plan  of  this  chapter  is  as  follows:  in  Section  2.1  we  summarize  both  the  theory 
of  Green’s  functions  for  first  order  linear  systems  and  the  theory  of  Newton  methods  for 
first  order  nonlinear  systems,  in  Section  2.2  we  develop  the  analytical  apparatus  to  be  used, 
and  in  Section  2.3  we  describe  the  numerical  schemes  themselves.  The  performance  of  the 
methods  is  illustrated  in  Section  2.4  with  numerical  examples. 

The  present  chapter  is  similar  to  [22]  in  that  while  it  is  based  on  a  sequence  of  fairly 
simple  observations,  the  details  of  the  algorithm  are  somewhat  involved.  We  attempt  in 
this  chapter  to  present  both  cursory,  qualitative  descriptions  as  well  as  detailed,  rigorous 
proofs. 


2.1  Mathematical  Preliminaries 

In  this  section,  we  summarize  the  relevant  properties  of  both  the  boundary  value  problems 
to  be  addressed  and  the  second  kind  integral  equations  to  be  used  for  their  solution.  Most 
of  the  results  are  classical  and  can  be  found,  for  example,  in  [11]  and  [13].  The  rest  are 
straightforward  generalizations  to  systems  of  ordinary  differential  equations  of  well-known 
facts  concerning  second  order  boundary  value  problems  (see,  for  example,  [14]). 

2.1.1  Notation  and  Definitions 


Definition  2.1  A  linear  first  order  system  of  ordinary  differential  equations  is  an  expres¬ 
sion  of  the  form 

$'(i) -f  p(i)  •  $(x)  = /(x),  (2.1) 

with  $  :  [a,c]  — ♦  R"  in  C^[a,c],  p  :  [a, c]  ^  L(R"’^”)  and  f  :  [a,  c]  — *  R"  continuous,  and 
L(Rnxn)  ^^fioting  the  linear  space  of  all  linear  operators  R”  — ►  R”. 

Definition  2.2  If  /(x)  =  0,  (2.1)  assumes  the  form 

$'(x)  +  p(x)  •  $(x)  =  0,  (2.2) 

and  is  referred  to  as  a  linear  homogeneous  first  order  system  of  ordinary  differential  equa¬ 
tions. 
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Definition  2.3  A  differentiable  function  $  :  [a,  c]  — >  R"  is  a  solution  to  a  linear  first  order 
boundary  value  problem  if  it  satisfies  an  equation  of  the  form  (2.1),  subject  to  boundary 
conditions  of  the  form 

A-$(a)  +  C-$(c)  =  7.  (2.3) 

with  6  L(R”^”),  and  7  €  R". 

Definition  2.4  If  ^  =  Q,  (2.3)  becomes 

>1  •  #(a)  +  C  ■  #(c)  =  0,  (2.4) 

and  is  referred  to  as  a  set  of  homogeneous  boundary  conditions. 

Definitions  2.5-2. 6  are  the  nonlinear  analogues  to  Definitions  2.1  and  2.3. 

Definition  2.5  A  nonlinear  first  order  system  is  defined  as  an  expression 

$'(!)  =  F($(x),x),  (2.5) 

with  $  :  [a,  c]  — ►  R”  in  C^[a,c],  F  :  R”'*'^  — >  R”  continuous. 

Definition  2.6  A  differentiable  function  $  :  [o,c]  — »  R”  is  a  solution  to  a  nonlinear  first 
order  boundary  value  problem  if  it  satisfies  an  equation  of  the  form  (2.5),  subject  to  boundary 
conditions  of  the  form 

>l-$(o)  +  C-$(c)  =  7,  (2.6) 

with  A,C  e  L(R"^"),  7  e  R”. 

Definition  2.7  A  continuous  function  G(x,t)  :  [a,c]  x  [a,c]  — >  L(R”^")  is  the  Green’s 
function  for  a  boundary  value  problem  (2.1),  (2.4)  if 

1.  is  continuous  except  at  x  =  t, 

2.  G{x  +  0,  x)  —  G{x  —  0,  x)  =  In  for  all  x  €  [a,  c], 

3.  ■^G{x,t)  +  p(x)  •  G{x,t)  =  0  for  all  x,t  €  [a,c],x  ^  t, 

4.  A  ■  G{a,t)  +  C  ■  G{c,t)  =  0  for  all  t  €  [a,  c]. 

Remark  2.1  Green’s  functions  are  the  principal  analytical  tools  which  enable  boundary 
value  problems  to  be  solved  as  second  kind  integral  equations.  However,  Green’s  functions 
are  known  or  computable  for  very  few  problems  (2.1),  (2.4).  Fortunately,  we  can  use  one 
of  the  known  Green’s  functions  when  constructing  the  second  kind  integral  equation  for  a 
particular  boundary  value  problem.  When  a  Green’s  function  unrelated  to  a  problem  (2.1), 
(2.4)  is  used  to  convert  that  problem  to  an  integral  equation,  we  will  refer  to  this  Green’s 
function  as  a  background  Green’s  function.  □ 
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Definition  2.8  A  function  T  :  [a,c]  — >  L(R"’‘")  is  called  a  fundamental  matrix  for  (2.2) 
if  it  is  nonsingular  and 

T'(*)  +  p(x)-T(i)  =  0  (2.7) 

for  all  X  €  [a,c]. 

We  define  boundary  condition  matrices  D  and  to  be  used  in  theorems  in  the  re¬ 
mainder  of  Section  2.1. 

Definition  2.9  Given  a  fundamental  solution  matrix  T  of  the  system  (2.2),  and  a  pair 
of  matrices  A,C  given  by  (2.3),  the  boundary  condition  matrices  D,Ds  G  L(R”’^")  are 
defined  by  the  formulae 


D  =  A  •  T(a) -1- C  •  T(c),  (2.8) 

=  A  +  C.  (2.9) 


We  define  a  residual  mapping  K  and  Newton  iterates  6k  to  be  used  in  a  Newton  method 
for  nonlinear  boundary  value  problems. 


Definition  2.10  Given  functions  Gq,Po  :  [a,c]  X  [a,c]  — >  L(R”^”),  we  define  the  residual 
mapping  K  :  R""*"^  — >  R"  by  the  formula 


K{(t{x),x) 


a{x)-po{x)-  f  Go{x,t)-cr{t) 

Ja 

-F  (  f  Go{x,i)  •  <T(t)dt,x') . 


dt 


(2.10) 


Definition  2.11  For  any  continuous  cq  •  [a,c]  -♦  R”,  U7e  refer  to  the  continuous  functions 
6k  :  [a,  c]  R”  as  Newton  iterates  if  for  each  k  =  1,2,..., 


6k  =  (Tk+i{x)  -  (Tk{x),  (2.11) 

with  each  continuous  ak  ■  [a,c]  — ►  R*^  recursively  defined  via  the  formula 

■  i(Tk+i{x)  ~  Ckix))  =  -K{(7kix),x),  k  =  0,l,---.  (2.12) 

Finally,  we  define  a  transposition  operator  to  be  used  in  the  chapter. 

Definition  2.12  Given  an  interval  [61,62]  C  R  and  an  operator  x  '•  ^'^[61,62]  ^  L(R"’^"), 
the  transpose  ■  (f'^[6i,  62])"  — >  R"  of  x  is  defined  by  the  formula 

X^i(^)=  f  Xit)-crit)dt,  (2.13) 

Jbi 

with  a  £  {L^Y . 
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2.1.2  Green’s  Functions  for  First  Order  Systems 

Theorems  2. 1-2.8  provide  the  tools  for  the  conversion  of  first  order  systems  of  differential 
equations  into  second  kind  integral  equations.  Theorems  2.1,  2.2,  2.5  and  2.6  are  well  known 
and  can  be  found,  for  example,  in  [11]  and  [13].  The  authors  failed  to  locate  the  remaining 
theorems  in  the  literature. 

Theorems  2. 1-2.2  provide  conditions  for  the  existence  and  uniqueness  of  solutions  to 
(2.1),  (2.4). 


Theorem  2.1  For  any  continuous  function  p  :  [a,  c]  — >  i(R"’‘"),  the  homogeneous  first 
order  system  (2.2)  has  exactly  n  linearly  independent  solutions. 

Theorem  2.2  If  the  matrix  D  defined  by  (2.8)  is  nonsingular,  then  there  is  a  unique  so¬ 
lution  $  to  the  equation  (2.1)  satisfying  homogeneous  boundary  conditions  (2.4).  Fur¬ 
thermore,  the  solution  to  the  homogeneous  equation  (2.2)  satisfying  homogeneous  boundary 
conditions  (2.4)  is  $(1)  =  0. 

The  purpose  of  the  following  two  theorems  is  to  permit  the  conversion  of  problems 
with  inhomogeneous  boundary  conditions  to  those  with  homogeneous  ones.  Theorem  2.3 
concerns  linear  problems  of  the  form  (2.1),  (2.3);  Theorem  2.4  concerns  nonlinear  problems 
of  the  form  (2.5),  (2.6). 

Theorem  2.3  If  the  boundary  condition  matrices  D,Ds  defined  by  (2.8),  (2.9)  are  both 
nonsingular,  then  the  solution  to  the  problem  (2.1),  (2.3)  is  given  by  the  formula 

$(x)  =  $(i)  +  I/,  (2.14) 

with  V  €  R”  given  by  the  formula 

i'  =  (^  +  C)-^-7,  (2.15) 

and  4  :  [a,c]  — >  R”  in  C^[a,c]  the  solution  to  the  first  order  system 

4'(i)  +  p(x)  •  #(ar)  =  f{x)  -  p{x)  ■  u,  (2.16) 

satisfying  homogeneous  boundary  conditions  (2.4)- 

Proof.  Since  the  matrix  D  is  nonsingular,  it  immediately  follows  from  Theorem  2.2  that 
there  exists  a  unique  $  satisfying  the  equation  (2.16).  Substituting  (2.14)  into  boundary 
conditions  (2.3),  we  obtain 


4  ■  ($(a)  +  I/)  +  C  •  ($(c)  +  i')  =  7.  (2-17) 

Now,  (2.15)  is  easily  obtained  from  the  combination  of  (2.17)  and  (2.4),  while  (2.16)  is  a 
result  of  substituting  (2.14)  into  (2.1). 
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Theorem  2.4  If  there  exists  a  solution  $  :  [a,c]  ^  R"  to  the  problem  (2.5),  (2.6),  and  if 
the  matrix  defined  by  (2.9)  is  nonsingular,  then  $  is  given  by  the  formula 


^(x)  =  $(x)  +  t/. 


with  1/  €  R"  given  by 


i/  =  (^+cr'.7, 

and  $  :  [a,c]  — >  R"  €  C^[a,c]  the  solution  to  the  nonlinear  boundary  value  problem 

=  F($  +  v,x) 


(2.18) 

(2.19) 

(2.20) 


with  homogeneous  boundary  conditions  (2.4). 


Proof.  Substituting  (2.18)  into  boundary  conditions  (2.6)  we  obtain 

A  •  ($(a)  +  I/)  +  C  •  ($(c)  +  I/)  =  7-  (2.21) 

Now,  (2.19)  is  easily  obtained  from  the  combination  of  (2.21)  and  (2.4),  while  (2.20)  is  a 
result  of  substituting  (2.18)  into  (2.5). 

Theorem  2.5  provides  an  explicit  construction  for  the  Green’s  function  for  a  boundary 
value  problem  with  a  known  fundamental  matrix  T.  Given  a  Green’s  function  for  a  homoge¬ 
neous  problem  (2.2),  (2.4),  Theorem  2.6  provides  an  explicit  solution  for  the  inhomogeneous 
problem  (2.1),  (2.4). 


Theorem  2.5  If  the  matrix  D  defined  by  (2.8)  is  nonsingular,  then  there  exists  a  unique 
Green’s  function  G  :  [a,c]  x  [a,c]  — ►  L(R"’^")  for  (2.2),  (2.4).  G  is  given  by  the  formula 


G{x,t) 


f  T(x).(T-Ht) +  </(«))  it<x), 
\T(x)-J(0  (t>x), 


(2.22) 


with  J  :  [a,c]  — >  L(R’*’'’*)  given  by  the  formula 


J{t)  =  -Z?-‘  •C-T(c)-T-n0- 


(2.23) 


and  T  ;  [a,c]  — »  L(R"’‘”)  the  fundamental  matrix  for  (2.2)  (see  Definition  2.8). 


Theorem  2.6  Given  a  Green’s  function  for  the  problem  (2.2),  (2.4),  the  solution  $  for  the 
problem  (2.1),  (2-4)  can  be  obtained  via  the  formula 

$(x)=  r  G{x,t)-  f{t)dt.  (2.24) 

Ja 

The  following  two  theorems  are  two  of  the  principal  analytical  tools  used  in  this  chapter. 
Theorem  2.7  is  used  to  reduce  a  linear  boundary  value  problem  (2.1),  (2.4)  to  a  second  kind 
integral  equation,  even  when  the  Green’s  function  for  the  problem  is  not  available;  Theorem 
2.8  is  used  in  the  same  fashion  to  reduce  nonlinear  boundary  value  problems  (2.5),  (2.4)  to 
nonlinear  second  kind  integral  equations. 
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Theorem  2.7  Suppose  po  :  [a,c]  — <■  L(R"’^”)  is  continuous,  To  :  [a,c]  — >  L(R"^”)  is  the 
fundamental  matrix  for  the  equation 

#'(a:)  +  po(x)  •  $(i)  =  0,  (2.25) 

and  Go  :  [a,c]  X  [a,c]  —>■  L(R"^")  is  the  Greenes  function  for  the  boundary  value  problem 
(2.25),  (2.4).  Suppose  further  that  the  matrix  D  defined  by  (2.8)  and  the  matrix  Dq  € 
L(Rnxn)  jgyjngj  }fy  (/jg  formula 

Do  =  A-  To(a)  +  C  •  To(c),  (2.26) 

are  both  nonsingular.  Then  the  solution  $  to  the  problem  (2.1),  (2.4)  can  be  obtained  via 
the  formula 

$(ar)=  r  Go{x,t)-a(t)dt,  (2.27) 

Ja 

with  (T  :  [a,  c]  — ►  R"  the  solution  to  the  second  kind  integral  equation 

a{x)  +  [p(i)  -  Pq{x)]  ■  f  Goix,  t)  ■  o{t)  dt  =  f{x).  (2.28) 

Ja 

Proof.  By  Theorem  2.2,  if  matrices  D,Do  are  nonsingular  then  the  problems  (2.1),  (2.4) 
and  (2.25),  (2.4)  have  unique  solutions,  and  therefore  the  background  Green’s  function  Go 
is  also  unique,  and  is  defined  by  Theorem  2.4.  Now,  (2.28)  is  obtained  by  substituting  (2.27) 
into  (2.1). 

Remark  2.2  If  po(x)  =  p(x),  then  the  solution  to  equation  (2.28)  is  trivially  <t  =  /. 
Our  working  assumption  is  that  for  some  background  problem  (2.25),  (2.4),  the  Green’s 
function  is  known  or  computable,  but  that  for  the  original  differential  equation  (2.1),  (2.4) 

the  Green’s  function  is  unavailable.  □ 

Theorem  2.8  Suppose  $  ;  [a,c]  — ►  R"  is  a  solution  to  (2.5),  (2.4)-  Suppose  further  that 
Po  :  [a,c]  — >  L(R"’'’*)  is  continuous,  and  To  ;  [a,c]  — >  L(R’*’'’*)  is  a  fundamental  matrix  for 
the  equation 

$'(i)  +  po(x)-$(x)  =  0,  (2.29) 

and  Go  '■  [a,c]  x  [a,c]  — ♦  L(R"’‘’*)  is  the  Green’s  function  for  the  boundary  value  problem 
(2.25),  (2.4).  Suppose  finally  that  the  matrix  Do  defined  by  the  formula 

Do  =  A-  To(a)  +  C  ■  To(c)  (2.30) 

is  nonsingular.  Then  $  can  be  obtained  via  the  formula 

$(x)=  r  Go{x,t)-cT{t)dt,  (2.31) 

Ja 

with  a  :  [a,c]  — ►  R”  the  solution  to  the  second  kind  integral  equation 

a(x)-po{x)-J  Goix,t)-<7{t)dt  =  f(^J  Goix,t)  ■  cr{t)dt,x j  .  (2.32) 

Proof.  Since  Do  is  nonsirgular,  the  background  Green’s  function  Go  is  unique,  and  there¬ 
fore  $  can  be  obtained  from  (2.31).  Now,  (2.32)  is  obtained  by  substituting  (2.31)  into 
(2.5). 
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2.1.3  Green’s  Functions  for  Particular  Equations 

Lemmas  2. 1-2.4  of  this  section  provide  fundamental  matrices  and  Green’s  functions  for 
two  particular  types  of  boundary  value  problems.  Lemmas  2.1,  2.2  are  easily  verified  by 
substituting  formulae  (2.34),  (2.35)  into  (2.7),  (2.22).  Similarly,  Lemmas  2.3, 2.4  are  verified 
by  substituting  formulae  (2.37),  (2.38)  into  (2.7),  (2.22). 

Lemma  2.1  A  fundamental  matrix  To  for  the  equation 


<&'(i)  =  0 


(2.33) 


is  given  by  the  formula 


To(x)  =  /„,  (2.34) 

with  n  the  dimensionality  of  the  problem  (2.33),  and  x  6  [a,c]  (in  accordance  with  standard 
practice,  /„  denotes  the  unity  operator  R**  —*■  R”^. 


Lemma  2.2  The  Green’s  function  Gq  corresponding  to  the  equation  (2.33)  subject  to  bound¬ 
ary  conditions  (2.4)  is  given  by  the  formula 


r  t.  /  4  -  (>1  +  C)-'  •  c 

Go(»,0=| 

{t  <  *), 

{t  >  x). 

(2.35) 

Lemma  2.3  For  any  X  E  R,  a  fundamental  matrix  Tq  for  the  equation 

$'(x)  +  A  •  $(x)  =  0 

(2.36) 

is  given  by  the  formula 

To(x)  =  •  In, 

(2.37) 

with  n  the  dimensionality  of  the  problem  (2.36),  and  x  G  [0, 1]. 


Lemma  2.4  The  Green’s  function  Gq  corresponding  to  the  equation  (2.36)  subject  to  bound¬ 
ary  conditions  (2.4)  is  given  by  the  formula 


Go{x,t) 


eMt-r)  .  _  eA(t-i) .  ^  g-A  .  Q-i 

-(A  +  e-^ -C)-! -C  (t>x). 


(2.38) 


2.1.4  Linear  IVansformations  for  Problems  with  Singular  Dq  or 

The  purpose  of  Theorem  2.9  is  to  permit  the  conversion  of  a  problem  (2.1),  (2.3)  to  a 
second  kind  integral  equation  (2.28).  For  most  problems,  Theorems  2.3  and  2.7  allow 
such  a  conversion,  but  Theorem  2.3  cannot  be  used  when  the  matrix  Ds  defined  by  (2.9) 
is  singular,  while  Theorem  2.7  cannot  be  used  when  the  matrix  Dq  defined  by  (2.26)  is 
singular.  We  remove  these  obstacles  in  this  section  by  providing  a  scheme  which  reduces  a 
problem  of  the  form  (2.1),  (2.3)  with  singular  matrices  Dq,Ds  to  a  problem  of  the  same 
form  with  nonsingular  Dq,  Ds- 

Theorem  2.10  generalizes  Theorem  2.9;  it  permits  the  conversion  of  nonlinear  problems 
of  the  form  (2.5),  (2.6)  to  nonlinear  integral  equations  of  the  form  (2.32). 
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Remark  2.3  If  only  the  matrix  Dq  is  singular,  one  can  always  choose  a  new  background 
Green’s  function  Go  for  which  Do  will  be  nonsingular.  However,  we  have  found  that  for 
most  problems  it  is  easier  to  develop  a  transformation  of  the  type  described  in  this  section 


than  to  develop  an  alternate  background  Green’s  function.  □ 

Theorem  2.9  Suppose  $  :  [a,c]  — >  R"  is  the  unique  solution  to  the  problem  (2.1),  (2.4). 
Suppose  further  that  To  :  [a,c]  — ♦  L(R"’^”)  is  a  fundamental  matrix  for  the  background 
equation  (2.25).  Suppose  finally  that  there  exists  $  :  [o,c]  —y  L(R"’'")  such  that  ^  £ 
G*[a,c],  det  0  for  all  x  6  [a,c],  and  the  matrix 

Do  =  A-  $(a)  •  To(o)  +  C  •  «(c)  •  To(c)  (2.39) 

is  nonsingular.  Then  the  equation 

r'(x)  +  •  (^'(x)  +  p(x)  .  $(x))  •  r(x)  =  9-\x)  ■  fix),  (2.40) 

subject  to  boundary  conditions 

A  ■  9(a)  .  r(a)  +  C  ■  9(c)  ■  r(c)  =  0.  (2.41) 

has  a  unique  solution  F  :  [a,c]  — ►  R",  and 

9(x)  =  ^(x)  •  r(x),  (2.42) 


for  all  X  £  [a,c]. 

Proof.  We  immediately  obtain  (2.41)  by  substituting  (2.42)  into  (2.4).  Now,  substituting 
(2.42)  and  its  derivative  into  (2.1),  we  get 

^'(x)  ■  r(x)  +  '^(x)  •  r'(x)  +  p(x)  •  «'(x)  •  r(x)  =  /(x),  (2.43) 

and  obtain  (2.40)  by  combining  (2.43)  with  the  fact  that  (x)  is  nonsingular  for  all  x  €  [a,  c]. 

Remark  2.4  Cle^ly,  the  transformed  problem  (2.40),  (2.41)  satisfies  the  conditions  of 
Theorem  2.7,  as  Do  defined  by  (2.39)  is  nonsingular.  However,  for  many  problems,  the 
boundary  condition  matrix  Df,  defined  by  the  formula 

Dy  =  A-  9(a)  +  C  •  $(c).  (2.44) 

is  singular,  and  therefore  the  transformed  problem  fails  to  satisfy  the  conditions  of  Theorem 
2  .t  If  one  needs  to  use  the  results  of  both  Theorem  2.7  and  2.3,  one  must  choose  a 
transformation  4’  such  that  both  Do  and  Df,  are  nonsingular. 

Of  course,  it  is  easier  to  choose  transformations  when  Do  =  I? at-  This  is  true  when  the 
background  Green’s  function  is  chosen  to  correspond  to  the  equation  9'  =  0.  By  Lemma  2.1, 
the  fundamental  matrix  for  this  equation  is  T  =  /„;  the  equivalence  for  this  fundamental 
matrix  of  (2.39)  and  (2.44)  is  readily  apparent.  □ 
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Theorem  2.10  is  the  nonlinear  analogue  of  Theorem  2.9;  the  proofs  of  the  two  theorems 
are  nearly  identical. 

Theorem  2.10  Suppose  $  :  [a,c]  — ►  R"  is  the  unique  solution  to  the  problem  (2.5),  (2.4). 
Suppose  further  that  To  :  [a,c]  — ♦  L(R"^”)  is  a  fundamental  matrix  for  the  background 
equation  (2.25).  Suppose  finally  that  there  exists  :  [a,c]  -»  L(R”^")  such  that  ^  G 
C^[a,c],  det  ^'(i)  ^  0  for  all  x  G  [a,c],  and  the  matrix 

Do  =  A-  ^{a)  ■  To(a)  +  C  ■  ^{c)  ■  To(c)  (2.45) 

is  nonsingular.  Then  the  equation 

r(i)  +  ■  r(i)  =  ^-\x)  ■  F{9{x)  ■  Tix),x)  (2.46) 

subject  to  boundary  conditions 

A  ■  ^(a)  -  r(a)  +  C  ■  ^(c)  •  r(c)  =  0.  (2.47) 

has  a  unique  solution  F  :  [a,c]  — ♦  R",  and 

^x)  =  9{x)-T{x),  (2.48) 


for  all  X  G  [<i,c]. 

2.1.5  Newton’s  Method  for  Nonlinear  Boundary  Value  Problems 

Theorems  2.4,  2.8  of  Section  2.1.2  reduce  nonlinear  boundary  value  problems  of  the  form 
(2.5),  (2.6)  to  nonlinear  second  kind  integral  equations  of  the  form  (2.32).  In  this  section, 
we  describe  the  convergence  properties  of  the  well-known  Newton’s  method  as  applied  to 
the  latter  (Theorem  2.12),  and  reduce  each  step  of  Newton’s  algorithm  to  the  solution  of  a 
linear  boundary  value  problem  of  the  form  (2.1),  (2.4)  (Theorem  2.11). 

Theorem  2.11  permits  each  Newton  iterate  Sk  defined  by  (2.11)  to  be  expressed  as  the 
solution  to  a  second  kind  integral  equation. 

Theorem  2.11  Suppose  K  :  R""''^  -♦  R”  in  C^[a,c\  defined  by  (2.10)  is  Frechet  differen¬ 
tiable  at  every  point  x  G  [a,c],  and  ^k  •  [a?*:]  — R”  is  defined  for  all  k  =  0, 1, . . .  via  the 
formula 

$t(i)  =  I  Go{x,t)-<Tkit)dt,  (2.49) 

Ja 

with  Ok  :  [a,c]  R"  defined  by  (2.12),  and  Gq  :  [a, c]  x  [o,c]  — ►  L(R"’^”)  the  Green’s 
function  for  (2.25),  (2.4).  Then  the  Newton  iterates  6k  '  [a,c]  — ♦  R"  G  C°[a,c]  given  by 
Definition  2.14  satisfy  the  equation 

6kix)  +  Dk(x)-  f  Go(x,t)-6k{t)dt  =  gk{x) 

J  a 


(2.50) 
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for  all  k  =  0, 1, . . .,  with  (1^  :  [a,  c]  — >  L(R"^”)  defined  for  all  k  =  0,1, .. .  by  the  formula 

^  ^  5F($*.(x),x)  ^  ^  /o  C1^ 

^ki=)  =  —  ~  ~ 

and  gk  :  [a,  c]  — >  L(R”^”)  defined  for  all  k  =  0,1,. . .  by  the  formula 

9kir)  =  Po(x)  •  ^a:)  +  F($it(x),  x)  -  (Tk(x).  (2.52) 

Proof.  (2.50)  is  obtained  by  substituting  the  Frechet  derivative  of  the  function  K  into 
(2.12),  and  substituting  (2.49),  (2.51),  (2.52)  into  the  resulting  equation. 

The  convergence  properties  of  Newton’s  method  have  been  thoroughly  studied.  Theorem 
2.12  is  one  fundamental  result,  and  can  be  found,  for  example,  in  [23]  (in  a  slightly  different 
form). 

Theorem  2.12  Suppose  $  is  the  unique  solution  to  (^-5),  (2.4),  ^  is  the  solution  to  (2.32), 
and  S  is  the  unique  solution  to  (2.50)  (so  that  the  linearization  (2.50)  to  the  equation  (2.32) 
is  nonsingular  at  a).  Then  there  exists  e  >  0  such  that  for  any  Oq  :  [a,c]  — >  R”  satisfying 
the  condition 

\\ao  -  ^11  <  e  (2.53) 

and  Newton  iterates  Ck  :  [a,c]  — ►  R"  defined  by  (2.11), 

1.  ||<7fc  -  d|l  <  e  for  all  k  =  1,2,..., 

2.  lim  Ok  =  d, 
fc-*oo 

3.  Ok  converges  to  a  quadratically. 

2.1.6  A  Lemma  from  Linear  Algebra 

Given  a  perturbation  of  the  unity  operator  I(i2)n  :  (T^)"  — »•  (T^)",  Lemma  2.5  provides  its 
inverse.  It  is  normally  used  when  the  rank  of  the  perturbation  is  low,  is  a  particular  case 
of  the  Sherman- Morrison  formula  (see,  for  example,  [19]),  and  is  easy  to  verify  directly. 

Lemma  2.5  For  any  two  vectors  U,V  ^  (L^)"’'"  such  that  -  U  ^  In, 

(/(i2)n  -  u  ■  +  U-{In-V^-  {/)-'  •  (2.54) 

2.2  The  Analytical  Apparatus 

In  the  remainder  of  this  chapter,  we  assume  that  the  solution  to  the  problem  (2.1),  (2.4) 
is  being  sought  on  the  interval  [a,c],  and  that  6  is  some  intermediate  point  (a  <  6  <  c). 
The  fundamental  observation  on  which  the  algorithms  of  Section  2.3  are  based  is  that  the 
solution  to  the  integral  equation  (2.28)  on  the  entire  domain  [a,  c]  can  easily  be  constructed 
from  the  solutions  of  two  independent  integral  equations,  one  defined  on  [a,  6]  and  one  on 
[6.c].  This  leads  naturally  to  a  recursive  algorithm,  in  which  independent  solutions  on  a 
large  number  of  subintervals  are  successively  merged  until  the  full  solution  is  obtained.  A 
precise  formulation  of  the  construction  and  the  resulting  numerical  scheme  will  require  some 
notation. 
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2.2.1  Notation 

We  will  denote  the  subintervals  [a,  6]  and  [6,  c]  of  [a,  c]  by  A  and  B,  respectively.  For 
convenience,  we  write  the  integral  equation  (2.28)  in  the  form 


a{x)-\-p{x)-  f  Go{x,t)-a{t)dt  =  f{x), 
Ja 


(2.55) 


with  p{x)  =  p{x)  —  Pq{x),  and  Gq  :  [a,c]  X  [a, c]  — *  L(R’*’'")  the  background  corresponding 
to  the  equation  (2.25)  subject  to  boundary  conditions  (2.4). 

We  define  the  operator  P  :  (Z/^[a,c])"  — ►  (L^[a, c])"  corresponding  to  (2.55)  by  the 
formula 

P{(r)(x)  =  <t(x) -{■  p{x)  ■  f  Go{x,t)‘Cr{t)dt,  (2.56) 

Ja 

so  that  we  have 

P(j  =  /.  (2.57) 

We  will  require  the  four  operators 

Paa  :  (Z^[a,5r-(l2[a,6])% 

Pab  :  {L^[b,c]r-^{L^[a,b]r  , 

Pba  :  (L^[a,b]r-*{L'^[b,c])\ 

Pbb  :  iLHb,c]r~^{L^b,c]r, 

defined  by  the  formulae 


PAA{<T)ix)  =  c{x)  +  pix)-  I  Goix,t)'<T{t)dt, 

Ja 

PabM{x]  =  p{x)-  j  GQ{x,t)-(T{t)dt, 

Jb 

PBA(<^)ix:)  =  P{x}-  f  Go{x,t)-(7{t)dt, 

Ja 


Pbb{(^){x) 


(7{x)  +  p{x)-J  Go{x,t)  •  (T{t)dt. 


(2.58) 

(2.59) 

(2.60) 
(2.61) 


We  define  the  operator  Q  :  (X^[a,  c])"’^”  (L^[a,c])"’^"  by  the  expression 


<3(x)(a:)  =  X(a:)  +  P{x)  •  f  Go{x,  t)  •  x(t)  dt. 

Ja 


(2.62) 


We  additionally  require  the  four  operators 


Qaa  :  (i^[a,6])”""-.(L2[a,6])"-", 
Qab  :  mb,c]r^^  ^  {L^[a,b]r^'^  , 
qba  : 

Qbb  :  {L^[b,c]r^^-^iL^{b,c]r^'^, 
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defined  by  the  formulae 

rh 

QAAix)ix)  =  Xix)  +  Pix)-  Goix,t)-x{t)dt,  (2.63) 

Ja 

Qab(x)(x)  =  t{x)-  j  Goix,t)-xit)dt,  (2.64) 

ri> 

Qba(x){x)  =  p{x)-  Go(x,t)-xit)dt,  (2.65) 

Ja 

Qbb{x){x)  =  x(x)  +  p{x)-  I  Go(x,t)-xit)dt.  (2.66) 

Jh 

We  also  require  the  functions  :  [a,c]  L(R”^'^)  defined  by  the  formulae 

^{x)  =  p(x)-To(a:),  (2.67) 

^i{i)  =  ^^0^(0  + -^0(0.  (2.68) 

u«(<)  =  Jo{t),  (2.69) 

with  To  the  fundamental  matrix  for  equation  (2.25),  and  Jo  :  [a,c]  L(R"^")  defined  by 

the  formula 

Jo(0  =  -V-C'-To(c).Toi(0,  (2.70) 

with  the  matrix  Dq  defined  by  (2.26),  and  the  matrix  C  given  by  (2.4). 


Given  a  function  /  6  (i^[a,c])’*,  we  will  follow  the  convention  of  denoting  its  restriction 
to  A  and  B  by  /|4  and  /|g,  respectively.  Similarly,  given  a  function  ip  €  (I^[a,  c])"^",  we 
will  denote  its  restriction  to  A  and  B  by  V’j/i  and  V’lB?  respectively.  Assuming  that  the 
operators  Paa,Pbb  are  nonsingular,  we  define  the  functions  77^  :  A  R",  tjb  :  B  — R” 
via  the  formulae 


=  PA\{f\A),  (2.71) 

nB  =  PbUAb)-  (2.72) 

Similarly,  assuming  that  the  operators  Q,QaAjQbb  are  nonsingular,  we  then  define  the 
operators 

X  :  [a,cr^"--L(R”=^"), 

<f>A  :  4”=^" L(R"’'"), 

4>b  :  -  L(R"’^”), 

via  the  formulae 

X  =  Q~\rP\  (2.73) 

=  QAAi^\A),  (2.74) 

=  QiBii’\B)-  (2.75) 
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Finally,  we  will  define  six  matrices  a^,a^,  af,af,  ai^,aR  €  L(R"’^")  by  the  formulae 


j 

/  Vi.{t)  ■  4>A{t)  dt  , 

'a 

(2.76) 

j 

f  i’r(0  -  <f>A{t)  di  , 

(2.77) 

af 

w 

1  M*)  ■  <t>B{t)  dt  , 

h 

(2.78) 

f  •  <^B(t)  dt  , 

'6 

(2.79) 

OtL 

= 

/  •  X(t)  dt  , 

'a 

(2.80) 

= 

f  Wr(0  •  X(0  dt  , 

'a 

(2.81) 

and  six  vectors  6i,S^,6^,S^,6i,6n 

eR 

via  the  formulae 

= 

/  VLit)-TjAit)dt, 

Ja 

(2.82) 

Si 

= 

f  Vftit)-T)A{t)dt, 

Ja 

(2.83) 

Sf 

= 

VL{t)-VB(t)dt, 

(2.84) 

S^ 

= 

VR{t)-r}Bit)dt, 

(2.85) 

Sl 

= 

/  Vi{t)-o{t)dt, 

Ja 

(2.86) 

= 

1  Vnit)  ■  ait)  dt, 

Ja 

(2.87) 

with  a  the  solution  to  equation  (2.57). 

2.2.2  Analysis  of  the  operators  Pab ^  Pba 

In  this  section,  we  observe  that  each  of  the  operators  Pab  and  Pba  is  of  rank 

n,  and  give 

simple  expressions  for  these  operators. 

Lemma  2.6  In  the  notation  of  the  preceding  section, 

Pab  =  t\A  ■  (2-88) 

Pba  =  i’lB  ■  (2-89) 

Proof.  VVe  obtain  (2.88)  by  observing  that  i  <  t  for  any  x  £  €  [6,c],  and  by 

using  (2. .59)  and  (2.69).  Similarly,  (2.89)  follows  from  the  combination  of  (2.60),  (2.69)  and 
observing  that  i  >  t  for  any  x  £  [6,c],  t  £  [a.fij. 
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Figure  2.1:  The  shaded  boxes  are  rank  n  matrices 

Remark  2.3  Just  as  we  decomposed  the  second  kind  operator  P  into  two  second  kind 
integral  operators  {PaAiPbb)  a^d  two  first  kind  integral  operators  {PaBiPba)i  can 
similarly  decompose  each  of  PaAiPbb  into  two  second  kind  integral  operators  and  two 
first  kind  integral  operators.  Clearly,  this  decomposition  can  be  applied  recursively  to 
yield  the  structure  shown  in  Figure  2.1,  in  which  each  block  diagonal  matrix  is  a  second 
kind  integral  operator,  while  the  remaining  block  matrices  correspond  to  rank  n  operators. 
This  immediately  suggests  a  recursive  algorithm  to  rapidly  obtain  solutions  to  the  integral 
equation  (2.57).  We  now  present  the  principal  lemmas  for  this  algorithm.  □ 


2.2.3  Recursive  solution  of  the  integral  equation 

We  now  consider  the  original  integral  equai'on  (2.57) 

Pa  =  f. 

The  main  result  of  this  section  is  the  following  lemma,  which  constructs  the  solution  a 
of  equation  (2.57)  from  of  equations  (2.71)  and  (2.72). 
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Lemma  2.7  If,  in  the  notation  of  Section  2.2.1,  all  six  operators  P,  Paa,  Pbb,Qi  Qaai 
Qbb  nonsingular,  and  the  matrices  Ai,  A2  €  L(R"’^")  defined  by  the  formulas 


(2.90) 

A2  =  /n-af-af,  (2.91) 

are  also  nonsingular,  then 

^lA  =  VA  +  d>A-  AJ'  •  (af  •  ),  (2.92) 

(Tib  =  T)s  +  d>B  ■  ■  (ot^  ■  -  S^).  (2.93) 

Proof.  Using  definitions  (2.56)  -  (2.61),  the  integral  equation 

Po  =  f  (2.94) 

can  be  rewritten  in  the  form 

P AAi<^\A)  +  PAB{cr\B)  =  /|^>  (2.95) 

Pba{<^\a)  +  ^bb(<7'ib)  =  f\B-  (2.96) 


The  expansions  (2.88)  and  (2.89)  for  Pab  and  Pba,  respectively,  can  then  be  used 
to  obtain  an  explicit  solution  to  the  coupled  equations  (2.95)  and  (2.96)  in  terms  of  the 
functions  r/^,  t]b,  (/>Ai  and  <f)B  defined  by  (2.71),  (2.72),  (2.74),  (2.75),  respectively.  Indeed, 
applying  the  operator  P^\  to  equation  (2.95)  and  the  operator  Pgg  to  equation  (2.96),  we 
have 


^\A  +  Paa  •  -P^b(<^|b)  =  PaHAa)^ 

(2.97) 

Pbb  ■  PBA{(r\A)  +  <^|B  =  -^BbI/ib)- 

(2.98) 

Substituting  (2.88)  and  (2.89)  into  (2.97)  and  (2.98)  yields  the  formulae 

^\A  +  Paa  •  '•^\A  '  ■  <^|B  =  VA, 

(2.99) 

Pbb  •  ^|B  •  vj  ■  (r\A  +  <^\B  =  VB, 

(2.100) 

or 

^\A  P  A  '  '  <^\B  —  ^Ai 

(2.101) 

T* 

4>B  -Vl  -(^lA  +<^\B  =  iB, 

(2.102) 

where  we  have  used  the  definitions  (2.74),  (2.75)  for  4>a  and  <t>Bi  respectively.  Now,  multi¬ 
plying  (2.102)  by  4>a  ■  t-’R  and  subtracting  it  from  (2.101),  we  obtain 


^A  '  ^’r  ■  4'B  ■  )  '  ^\A  —  V.A  4*A  '  ^R  ■  VB- 


(2.103) 
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Similarly,  multiplying  (2.101)  by  4>b  •  subtracting  it  from  (2.102)  results  in  the 


equation 

(/(£,2)n  -  <j>B  ■  vj  ■  <t>A  ■  vl)  ■  £T|b  =  nB  -  <t>B  •  •  VA-  (2.104) 

Due  to  (2.76),  (2.79),  and  (2.82),  (2.85)  we  can  rewrite  these  equations  in  the  form 

(/(i2)n  -  <t>A  ■  (of  •  uj))  •  =  Va-4>a-S^,  (2.105) 

(^(L2)"  -  4>b  ■  («f  •  »D)  •  o’lB  =  Vb-4>b  (2.106) 

By  application  of  Lemma  2.5,  we  obtain 

(^\A  =  4-  <f>A  •  {In  -Otf-vl-  (f>A)~^  •  of  •  vj)  •  {tJa  -  <f>A  '  ),  (2.107) 

^|B  -  +  4>B  ■  (In  -O-i  4>b)~^  •  O'i  •  vj)  •  {t)b  -  4>B  ■  )•  (2.108) 


The  equations  (2.92),  (2.93)  are  now  obtained  from  equations  (2.107),  (2.108)  and  equations 
(2.90),  (2.91). 

Remark  2.6  Suppose  that  6i  and  62  aje  a  pair  of  real  numbers  such  that  a  <  bi  <  b2  <  c, 
and  that  the  interval  [61,62]  is  denoted  by  C.  We  will  denote  by  Pcc  the  restriction 
of  the  operator  P  to  the  interval  C,  and  denote  by  Qcc  the  restriction  of  the  operator 
Q  to  the  interval  C.  Assuming  that  PcciQcc  nonsingular,  we  define  the  functions 
7JC  :C  -*R^,4>c  :C  L(R"^”)  by 

ric  =  PchiAc),  (2.109) 

<f>C  =  Qd'Mc)-  (2.110) 

By  applying  the  above  lemma  twice  (once  for  the  subinterval  [a,  61]  and  once  for  [0,62]),  we 
may  ea.sily  observe  that  there  exists  A  6  R”  such  that 

cr(a:)  =  7/c(ar)  +  (?!>c(x)  •  A  (2.111) 

for  all  X  €  C.  The  exact  expression  for  the  vector  A  is  complicated,  but  irrelevant  for  the 
purposes  of  this  chapter.  The  existence  of  a  relation  of  the  form  (2.111),  however,  wifi  be 

critically  important  in  Section  2.3.  □ 

The  following  corollary  constructs  the  solution  x  of  equation  (2.73)  from  4>a,4>b  of 
equations  (2.74)  and  (2.75). 

Corollary  2.1  If,  in  the  notation  of  Section  2.2.1,  all  six  operators  P,  Paa>  PbBi  Q>  Qaa, 
Qbb  are  nonsingular,  then 

X|4  =  <^^-Aj'.(/„-of),  (2.112) 

X|B  =  <Ab  •  Ar' •  (/„  -  Q^).  (2.113) 

with  the  matrices  af  and  af  defined  by  equations  (2.76)  and  (2.79),  and  the  matrices 
Ai,A2  defined  by  equations  (2.90)  and  (2.91). 
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Proof.  Substituting  in  equations  (2.107),  (2.108)  the  functions  <f>A,<pB  defined  by  (2.74), 
(2.75)  for  the  functions  rjA,T)B  defined  by  (2.71),  (2.72),  and  the  matrices  defined 

by  (2.76),  (2.79)  for  the  vectors  defined  by  (2.82),  (2.85),  we  obtain 

X|.4  =  4>a-<I>a-ccr+<I>a-^2^ -af-ai -<f>A-^2^ -af,  (2.114) 

X|B  =  4>b  -  4>b  ■  Oif  +  (f>B  ■  •  oil  ■  otf  -  4>b  ■  •  Oti  ■  af  •  Ol  ■  (2.115) 

The  expressions  (2.112),  (2.113)  are  now  easily  obtained  from  the  equations  (2.114),  (2.115). 

2.2.4  Further  Analytical  Results 

We  now  collect  a  number  of  identities  which  are  necessary  for  Algorithm  A,  to  be  presented 
in  Section  2.3. 

Corollary  2.2  provides  analytical  expressions  for  the  inner  products  and  6r  defined 
by  (2.86),  (2.87)  in  terms  of  the  restricted  inner  products  and  6^  defined  by 

(2.82)-(2.85). 

Corollary  2.2  If,  in  the  notation  of  Section  2.2.1,  all  six  operators  P,  PaAi  Pbb-,  Qi  Qaa, 
Qbb  nonsingular,  then 

=  I  Vi(t)  •o(t)dt  =  [  Ui(t)-0|^(t)dt+  f  Vi,{t)-CT\Bit)dt 

Ja  Ja  Jh 

=  +  +  (2.116) 
=  y  Vfiit)  ■  a(t)  dt  =  J  VR{t)-a\A{t)dtA  !;«(<)•  cr|B(t)  dt 

=  +  +  a^A2‘•(af +  •A^'•(a^^f-^^).  (2.117) 

Proof.  Multiplying  equation  (2.92)  by  vj  and  uj,  and  equation  (2.93)  by  vj  and  we 
obtain 

f\,{t)-o^A{t)dt  =  (2.118) 

Ja 

rv,it)-a^B{t)dt  =  +  af  •A^‘•(a^^f-^^),  (2.119) 

J  b 

l\Rit)-o^^{t)dt  =  S^  +  ai-A^^-iaf-S^-Sf),  (2.120) 

J  a 

J\Rit)-a^Bit)dt  =  +  af  •A^'•(a^^f -5^).  (2.121) 

Now,  expressions  (2.116),  (2.117)  are  easily  obtained  from  (2.118)-(2.121). 

CoroUary  2.3  is  similar  to  Corollary  2.2,  but  uses  Xi  tbe  matrix  valued  function  defined 
by  (2.73),  in  place  of  a,  the  vector  valued  function  defined  by  (2.57).  While  the  two 
corollaries  concern  different  objects  (the  vectors  6i,,6r  in  CoroUary  2.2,  the  matrices  ai,QR 
in  CoroUary  2.3),  their  proofs  are  nearly  identical. 
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Corollary  2.3  If,  in  the  notation  of  Section  2.2.1,  all  six  operators  P,  Paa,  Pbb,  Q,  Qaa, 
Qbb  ore  nonsingular,  then 

(^L  =  j  V:,{t)-x(t)dt  =  j  v^it)  ■  x\Ait)  dt  +  Vi,it)-X\Bii)dt 

=  of  •  A2 '•(/„- of) +  af-Ar'-(/„-af),  (2.122) 

Or  =  f  Vii(t)  •  x(0  dt=  f  VR(t)  -  X|yi(0  dt+  f  VR{t)  ■  X|s(0  dt 
J  a  J  d  Jh 

=  ai  •  AJ'  •  (/„  -  of)  +  of  •  Afi  •  (/„  -  a^).  (2.123) 

Finally,  combining  Lemma  2.7  with  the  expressions  (2.112)-(2.113),  we  have 

Corollary  2.4  Suppose  that  in  the  notation  of  Section  2.2.1,  all  six  operators  P,  Paa, 
Pbb,  QiQaa,Qbb  nonsingular.  Suppose  further  that  the  function  F  is  defined  by  the 
formula 

F(x)  =  x-X  +  CT.  (2.124) 

with  X  €  R".  Then  on  the  interval  [a,<>], 

F{x)  =  d>A{x)  •  p  +  gA{x),  (2.125) 

with  p  6  R"  defined  by  the  formula 

p  =  Aj ‘(A  -  of  •  (A  -  -  6f).  (2.126) 

Similarly,  on  the  interval  [6,  c], 

Fix)  =  <t>B{x)  ■  +  tib{x),  (2.127) 

with  V  €  R"  defined  via  the  formula 

1/  =  AJ-*  •  (A  -  a^  (A  -  6f)  -  6t).  (2.128) 

Proof.  Restricting  (2.124)  on  the  subintervals  A,  B  of  [a,c],  respectively,  we  have 

F\a  =  Xiyi'A  +  a|^,  (2.129) 

F\b  =  X\b  '  A  +  <xiB-  (2.130) 

Combining  (2.129),  (2.130)  with  (2.92),  (2.93),  (2.112),  (2.113),  we  obtain 

F^a  =  </>^-Aj'-(/„-af).A  +  (,74  +  <^4-A2'-(Qf -^f-^f)),  (2.131) 

F\b  =  ^B•Al-'•(/„-af).A  +  (77H  +  <^B.A^'.(a^^f-^^)).  (2.132) 

Now,  the  expressions  (2.126),  (2.128)  immediately  follow  from  the  comparisons  of  (2.125), 
(2.127)  with  (2.131),  (2.132),  respectively. 
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2.3  Description  of  the  Algorithms 

We  turn  now  to  the  construction  of  the  fast  algorithm  for  the  solution  of  the  integral 
equation  (2.57) 

based  on  the  apparatus  developed  in  Section  2.2.  The  main  tool  at  our  disposal  is  the  ability 
to  merge  the  solutions  of  restricted  versions  of  the  integral  equation  in  adjacent  subintervals 
(Lemma  2.7).  As  this  suggests  a  recursive  procedure,  we  begin  by  subdividing  the  whole 
interval  [a,  c],  on  which  the  solution  to  (2.57)  is  sought,  into  a  large  number  of  subintervals. 
For  the  sake  of  simplicity,  we  assume  that  m  is  a  positive  integer  and  that  Af  =  2”*  is  the 
number  of  subintervals  created.  The  boundary  points  of  the  subintervals  are  then  defined 
by  a  strictly  increasing  sequence  of  numbers 

61,62,  ...,6m,  ^M+ij  (2.133) 

with  61  =  a  and  6m+i  =  c.  For  each  i  =  we  define  the  interval  B^  via  the 

expression 

Br  =  [h.,h.+i],  (2.134) 

and  create  a  hierarchy  of  intervals  by  recursively  merging  adjacent  pairs.  That  is,  for 
each  j  =  m  —  1, . . .,  1,0,  and  t  =  1, . . M,  we  define 

=  (2-135) 

We  will  refer  to  each  fixed  /  as  a  level.  We  will  also  refer  to  the  two  intervals  B^tii  and 
as  children  and  to  the  larger  interval  5,-  as  a  parent. 

It  is  obvious  that 

Bi  =  [5i4.(i_i).2">-<  -.217.-1],  (2.136) 

and  that  for  each  level  /, 

2' 

[a,c]=UBi.  (2.137) 

:=1 

2.3.1  Notation 

Generalizing  the  notation  of  Section  2.2,  we  will  denote  by  Pij  the  restriction  to  the  interval 
BI  of  the  integral  operator  P,  so  that 

Pij{a){x)  =  (t{x)  +  P(x)-  Go{x,t)  ■  cr{t)dt  (2.138) 

2”>-< 

for  any  a  €  L^{B[)^.  Similarly,  we  will  denote  by  Qi^i  the  restriction  to  the  interval  5-  of 
the  integral  operator  Q,  so  that 

=  Xix)  + I  GQ{x,t)-x{t)dt 


(2.139) 
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for  any  x  €  For  each  £,•  we  wUl  define  the  functions  t/,,/  :  5,-  ^  R” ,  :  R,-  — ' 

L(Rnxn)  3^  solutions  of  the  equations 

P,.i{Vi,i)  =  /|B«.  (2-140) 

(2-141) 

provided  th^  operators  Pi^t,Qij  are  nonsingular. 

Remark  2.7  Suppose  now  that  the  operators  Pij,Qij  are  nonsingular  on  the  interval  R-. 
Then,  due  to  (2.111),  there  exists  A‘’^  €  R"  such  that 

ct{x)  =  +  4>M  ■  (2-142) 

for  all  I  €  R-.  D 


For  each  /  =  0, 1, . . .,  m,  and  i  =  1,2, ..  .,2*,  we  define  the  matrices  a'l  ,0'^  E  L(R”’‘”) 
by  the  formulae 


and  the  vectors 


i.l 

7*’>  +  .2'"-' 

ij 

/S+»-2’”-- 

S'A  €  R" 

by 

the  formulae 

jc‘4 

+  .  2’"-' 

Vi 

-'*l+(i-l)  2">-< 

Sr 

= 

+  .  2”*-' 

2'"-' 

(2.143) 

(2.144) 


(2.145) 

(2.146) 


2.3.2  Discretization  of  the  Restricted  Integral  Equations 


Choosing  an  integer  p  >  1,  we  construct  the  p  scaled  Chebyshev  nodes 


(2j  -  l)x 

2p 


+ 


^t+i  4- 


j  =  l,2,...,p  (2.147) 


on  each  of  the  intervals  i  =  1,2, ...,M.  We  then  discretize  the  two  integral  equations 
(2.140),  (2.141)  via  a  Nystrom  algorithm  based  on  p-point  Chebyshev  quadrature  (see,  for 
example,  [22]).  The  resulting  approximations  to  the  functions  77,,/,  4>i^i  at  the  nodes  t/  will 
be  denoted  by 

Vx.i  =  , 

respectively. 
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Remark  2.8  It  is  well-known  that  the  order  of  convergence  of  the  approximations 
to  the  functions  t?,,/,  <f>ij  is  p.  Since  all  subsequent  steps  in  the  construction  of  an  approximate 
solution  a  to  the  integral  equation  (2.57)  are  analytic,  the  convergence  rate  of  the  full 
algorithm  depends  entirely  on  the  parameter  p.  For  example,  by  using  16  scaled  Chebyshev 
points  on  each  subinterval  at  the  finest  level,  one  obtains  a  sixteenth  order  method. 

For  many  boundary  value  problems  (2.1),  (2.3),  the  matrix  function  p  :  [a,c]  — ►  L(R”^”) 
is  singular  in  the  interval  [a,c]  (see  Example  2.3  in  Section  2.4  below).  For  such  problems, 
Chebyshev  quadrature  will  yield  acceptable  convergence  only  if  the  singularity  is  removed 
by  means  of  an  appropriate  choice  of  the  “background”  Green’s  function.  □ 

Remark  2.9  The  algorithm  of  this  section  makes  extensive  use  of  the  apparatus  of  Cheby¬ 
shev  interpolation,  quadratures,  composite  quadratures,  etc.  This  apparatus  is  quite  well- 
developed,  and  can  be  found  in  various  forms  in  [16],  [18],  [20].  For  a  detailed  description 
in  the  form  most  convenient  for  our  purposes,  we  refer  the  reader  to  [22].  □ 


2.3.3  Informal  Description  of  the  Algorithm  for  Linear  DDEs 

We  begin  by  directly  solving  the  two  integral  equations  (2.140),  (2.141)  on  each  subinterval 
BJ"  at  the  finest  level,  as  discussed  in  the  preceding  section.  Equation  (2.142)  then  shows 
that  <7  restricted  to  B-"  can  be  expressed  as  a  linear  combination  of  the  two  solutions 
Thus,  it  remains  only  to  determine  the  coefficients  A*’”*  6  R”  for  each  of  the 
A/  subintervals  i?"*.  Fortunately,  this  can  be  done  recursively.  To  see  this,  suppose  that, 
at  some  coarse  level  /  <  m  -  1,  we  are  given  the  coefficient  A’’^  for  the  subinterval  Bj. 
Then  Corollary  2.4  provides  formulae  for  the  calculation  of  the  corresponding  coefficients 

Jqj.  intervals  and  B^f^,  respectively.  On  the 

coarsest  level,  we  observe  that  A°’*  =  0,  i.e.  the  solution  of  equation  (2.140)  on  the  whole 
interval  [a,c]  is  simply  ct. 

However,  the  formulae  (2.126)  and  (2.128)  of  Corollary  2.4  contain  the  matrices  a^' 

and  the  vectors  These  quantities 

are  also  computed  recursively  but  in  the  opposite  direction,  namely,  from  the  finest  level  to 
the  coarsest.  They  are  certainly  available  at  level  m  directly  from  the  definitions  (2.143)- 
(2.146).  For  the  interval  B\  at  any  coarser  level  /  <  m  -  1,  Corollaries  2.7  and  2.1  describe 
how  matrices  a'l  and  vectors  are  obtained  from  the  matrices  and  vectors 

Si^Sp,  of  the  two  child  intervals. 

To  summarize,  the  algorithm  consists  of  three  parts.  First,  a  sufficiently  fine  subdivision 
61,62. ..  .,5.vf-n  of  the  inter\'al  [a,c]  is  chosen  so  that,  on  each  of  the  intervals  the 

functions  can  be  accurately  represented  by  a  low  order  Chebyshev  expansion. 

On  each  of  the  intervals  B,_m,  the  equations  (2.140)-(2.141)  are  solved  (approximately)  by 
direct  inversion  of  the  linear  system  arising  from  a  Nystrom  discretization.  Second,  the 
matrices  o'/.Qr^  and  vectors  are  computed  in  an  upward  sweep,  beginning  at  the 

finest  level  m.  Finally,  the  coefficients  A'-^  are  computed  in  a  downward  sweep,  beginning 
at  the  coarsest  level.  The  desired  function  a  is  then  recovered  on  each  subinterval  from 
equation  (2.142). 
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The  following  is  a  more  detailed  description  of  the  numerical  procedure. 

Algorithm  A 

Comment  [Define  the  computational  grid.] 

Create  M  =  2"'  subintervals  on  [a,  c]  by  choosing  a  sequence  of  boundary  points  6i,  62, . . . , 
f>M,f>M+i  with  bi  =  a  and  6^+1  =  c.  Choose  the  number  p  of  Chebyshev  nodes  on  each 
interval  B["  =  [6i,  hi+i]  for  i  =  1, . . . ,  M.  Determine  the  locations  of  the  scaled  Chebyshev 
nodes  r/,  ,  rf  on  each  interval  B[",  and  evaluate  the  functions  /,  ip  at  these  nodes, 

obtaining 

Step  1. 

Comment  [Construct  the  approximate  solutions  on  each  interval  BJ”.] 

do  i  =  1,2, . . M 

(1)  Construct  the  two  p  ■  n  y.  p  ■  n  linear  systems  on  Bj  obtained  through  a  Nystrom 
discretization  of  the  corresponding  integral  equation. 

(2)  Solve  the  two  p  ■  n  y  p  ■  n  linear  systems  on  B[  by  Gaussian  elimination, 
obtaining  the  values 

end  do 


Step  2. 

Comment  [Construct  the  matrices  Ot'" ,  Or’"  and  vectors  6],''" ,  6'/^  on  each  interval  B["  at  the 
finest  level.] 

do  i  =  1,2, . . .,  M 

Evaluate  the  matrices  and  vectors  51’’",  using  the  p— point 

Chebyshev  quadrature  formula, 
end  do 


Step  3  (Upward  Sweep). 

Comment  [Construct  the  matrices  oI'^Or*  and  vectors  51’*,  5r*  for  all  intervals  at  all  coarser  levels 

/  =  m  —  l,m  -  2, . .  .,0  ] 


do  1=  m-1,  0,  -1 
do  i=l,  2* 

Compute  the  matrices  q1’*  ,  Or*  and  vectors  51’* ,  5r*  from  the  corresponding  data  in 


end  do 
end  do 


the  two  child  intervals  (af 


2t-l.J+l  ^2i-l,J+l  ^2i,l+l  g2i-l,/+l  _(;2i-l,l+l 


Sr  ^  ),  using  the  results  of  Corollaries  2.2  and  2.3. 
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Step  4  (Downward  Sweep). 

Comment  [Construct  the  coefficient  A*’"*  for  all  intervads  at  the  finest  level.] 

Set  =  0. 

do  l=0,m-l 
do  i=l,  2^ 

Use  Corollary  2.4  to  compute  the  coefficients  A*"*"*’^  A^*  *"''^, 

for  the  child  intervals  and  from  the  coefficient  A’’*  of  the  parent  interval 

end  do 
end  do 


Step  5. 

Comment  [Compute  the  solution  cr  of  equation  (2.57)  at  the  nodes  . .  .,Tf  for  each  interval 

S'"  at  the  finest  level.] 

do  i=l,  M 
doj=l,p 

Determine  the  values  of  the  solution  a  of  equation  (2.57)  at  the  node  rf  via 
formula  (2.142). 
end  do 
end  do 


Step  6. 

Comment  [Compute  the  solution  0  of  equation  (2.1)  from  the  values  of  a] 

Evaluate  the  integral  (2.27),  by  using  composite  Chebyshev  quadrature 
(see  Remark  2.11  below). 

Remark  2.10  Inspection  of  the  above  algorithm  shows  that  the  amount  of  work  required 
is  of  the  order  0{M  •  •  n^).  Step  1  involves  solving  two  (p  x  n)  x  (p  x  n)  linear  systems  for 

each  of  the  M  intervals.  Steps  2-5  require  no  more  than  0{M -p-n^ -{log  p-\-n))  operations. 
Since  N  =  M  •  p  is  the  total  number  of  nodes  in  the  discretization  of  the  interval  [a,c], 
we  can  write  the  CPU  time  estimate  in  the  form  0{N  ■  p^  •  n^).  The  cost  of  evaluating 
the  solution  $  of  the  differential  equation  (2.55)  from  the  integral  representation  (2.27)  is 
0{N  ■  logp  •  n)  (see  Remark  2.11  below).  □ 

Remark  2.11  The  final  step  in  the  algorithm  involves  the  evaluation  of  an  integral  of  the 
form  (2.27)  at  each  of  the  Chebyshev  nodes  on  each  subinterval  B^',  namely 

=  /  Go(Ti,t)-(T{t)dt. 

Ja 


(2.148) 
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If  these  integreds  were  calculated  independently  for  each  rf ,  the  amount  of  work  required 
would  be  of  the  form  0(N^  •  n),  and  would  dominate  the  construction  of  the  function  In 
fact,  this  is  unnecessary,  for  we  may  write 

^(rf)  =  T(Tf)-  r  (2.149) 

Ja  Jbi 

+  /  Ur(0  •  (^(t)  dt+  f  VR{t)  ■  a{t)  di  , 

where  we  have  used  the  representation  (2.22)  and  the  fact  that  rf  lies  in  the  interval  = 
Step  6  can  then  be  written  in  detail  as  follows: 

Step  6  (a). 

Comment  [Precompute  the  integrals  of  Vl  ■  a  and  Vr  •  <t  on  each  subinterval  SJ"  by  Chebyshev 
quadrature.  These  integrals  will  be  denoted  and  /«,  respectively.) 

do  i=l,  M 

=  v,{t)<r(t)dt. 

wr)  =  ft‘  v„(t)-a(t)dt. 

end  do 

Step  6  (b). 

Comment  [March  across  interval  from  a  to  c,  computing  #  at  each  node  in  the  discretization.  The 
variables  Jt,  and  will  be  used  to  eiccumulate  the  integrals  '  Vc(t)  ■  (r(<)  dt  and  v„(t)  ■  a(t)  dt, 
respectively. 

Set  Jr  =  I,(Br). 

Set  =  0. 

do  i=l,  M 
doj=l,p 

For  each  rf ,  compute 

=  r(i )  [ j, + /;/  v,(t)  ■  <T(t) dt + v„(t)  ■  a(t)  dt + j„] 

end  do 

Jl  =  Jl  +  hiBD 
Jr  =  J«  -  /R(Sn.i) 
end  do 

Thus,  the  amount  of  work  required  in  Step  6(a)  is  0{N  ■  n).  The  integrals  required  on 
each  subinterval  in  Step  6  (b)  can  be  computed  by  spectral  integration  (see,  for  example, 
[22])  using  0{p  ■  logp  ■  n)  work.  The  total  cost  is  therefore  of  the  order  0{M  •  p  ■  logp  •  n) 
or  0{N  ■  logp  •  n).  □ 


28 


CHAPTER  2.  TWO-POINT  BOUNDARY  VALUE  PROBLEMS 


2.3.4  Informal  Description  of  a  Simplified  Algorithm  for  Linear  ODEs 

Figure  2.1  suggests  a  simpler  version  of  Algorithm  A,  in  which  the  solutions  are 

obtained  for  all  levels  I  =  m,Tn—  1,...,0,  with  the  solution  a  =  7/0,1.  The  time  complexity 
of  the  resulting  algorithm  is  0{N -q^ -k- N  '\og.^{N)-(^)\  since  this  compares  unfavorably 
with  Algorithm  A,  we  merely  sketch  the  simplified  algorithm  below. 

Algorithm  A' 

Comment  [Define  the  computational  grid.] 

Create  M  =  2”*  subintervals  on  [a,  c]  by  choosing  a  sequence  of  boundary  points  61, 62, ... , 
with  61  =  a  and  buf+i  =  c.  Choose  the  number  p  of  Chebyshev  nodes  on  each 
interval  =  [6i,  for  i  =  1, . . . ,  M.  Determine  the  locations  of  the  scaled  Chebyshev 
nodes  rl,  r?, . . . ,  rf  on  each  interval  S[",  and  evaluate  the  functions  /,  16  at  these  nodes, 
obtaining 


Step  1. 

Comment  [Construct  the  approximate  solutions  on  each  interval  B-”.] 

do  i  =  1,2,  ...,A/ 

(1)  Construct  the  two  p  -  n  x  p  n  linear  systems  on  B'  obtained  through  a  Nystrom 
discretization  of  the  corresponding  integral  equation. 

(2)  Solve  the  two  pnx  pn  linear  systems  on  B\  by  Gaussian  elimination, 
obtaining  the  values 

end  do 


Step  2. 

Comment  [Construct  the  matrices  and  vectors  on  each  interval  at  the 

finest  level.] 

do  i  =  1, 2, . . .,  M 

Evaluate  the  matrices  and  vectors  using  the  p— point 

Chebyshev  quadrature  formula, 
end  do 


Step  3  (Upward  Sweep). 


Comment  [Construct  the  solutions  T]i,i,<i>i,i,  matrices  and  vectors  for  all  intervals 

at  all  coarser  levels  /  =  m  —  1,  tti  —  2, . . . ,  0.] 


do  1=  m-1,  0,  -1 


do  i=l,  2 


1 


(1)  Compute  the  matrices  and  vectors  from  the  corresponding  data 


2i-l,/+l  2i-l,J+l  ^2«.(+l 


2»,(+l  r2«-l,(+l  f2i-l,/+l 
1 


in  the  two  child  intervals  (q 

using  the  results  of  Corollaries  3.4  and  3.5.  Alternatively,  compute 
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the  matrices  a]’' vectors  via  formulae  (2.143)-(2.146)  using  the  p-point 

Chebyshev  quadrature  formula. 

(2)  Compute  the  solutions  from  data  in  the  two  child  interveils  (>72«_i,/+i! 

using  the  results  of 

Lemma  2.7  and  Corollary  2.1. 
end  do 
end  do 


Step  4. 

Comment  [Compute  the  solution  <j>  of  equation  (2.1)  from  the  values  of  <t.] 
Evaluate  the  integral  (2.27),  using  composite  Chebyshev  quadrature. 


2.3.5  Informal  Description  of  the  Algorithm  for  Nonlinear  DDEs 

The  nonlinear  algorithm  is  a  straightforward  application  of  Algorithm  A  described  in  Section 
2.3.3.  The  solution  is  obtained  using  Newton’s  method  for  nonlinear  ODEs;  each  Newton 
iterate  is  obtained  by  solving  the  linearized  problem  (2.50)  via  Algorithm  A. 

As  with  Algorithm  A,  we  subdivide  the  interval  [a,  c]  into  a  large  number  of  subintervals 
M;  for  simplicity  we  assume  M  =  2’",  with  m  a  positive  integer.  As  before  the  boundary 
points  bi,b2,...,bM,bM+\  are  defined  by  (2.133),  and  the  intervals  Bj,(l  <  I  <  m),(l  < 
i  <  2')  by  (2.134). 

On  the  step  of  the  Newton  process,  Algorithm  A  is  applied  to  the  integral  equation 

P%  =  9k,  (2.150) 

with  the  operator  :  (T^[a,c])"  — +  (T^[a.c])"  defined  by  the  formula 

P*(4)(z)  =  h{x)  +  fifc(ar)  •  r  Go{x,  t)  ■  6k{t)  dt,  (2.151) 

Ja 

with  6k  the  solution  of  the  integral  equation  (2.50),  ftfc  given  by  (2.51),  and  gk  given  by 
(2.52).  The  integral  equations  (2.140),  (2.141)  now  assume  the  form 

PtiiVij)  =  gkiB>.  (2-152) 

(2.153) 

with  the  operator  P-*',  :  (T^[a,c])"  — ‘  (Z/^[a,c])’‘  defined  by  the  formula 

Pl^_,i6k)ix)  =  6k{x)  +  Qk{x)-  Go{x,t)-6kit)dL  (2.154) 

and  the  operator  Q^i  :  (Z/^[a,  c])"’^"  — "  (Z/^[a,c])’*’'’'  defined  via  the  formula 

Qu(\)(x)  =  X(^)  +  Qkix)  ■  G'o(x,  t)  ■  x(t)  dt.  (2.155) 
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Once  Algorithm  A  has  computed  the  solution  Sk  to  (2.150),  we  obtain  Uk+i  via  (2.11),  and 
via  (2.49). 

The  nonlinear  algorithm  requires  an  initial  approximation  $0^0  solution  $  and 

its  derivative  of  equation  (2.5),  and  we  assume  that  both  are  supplied  by  the  calling 
program,  oq  is  obtained  from  $0,^0  the  identity 

CTo(x)  =  $o(x)  +  po(a:)  •  ^x). 

The  procedure  is  terminated  when  the  stopping  criterion 

<  e  (2.157) 

is  satisfied,  with  e  provided  by  the  calling  program.  Since  Newton’s  method  frequently  fails 
to  converge,  the  calling  program  also  permits  a  certain  majdmum  number  of  iterations,  after 
which  the  algorithm  stops,  signaling  failure. 

The  following  is  a  more  detailed  description  of  the  numerical  procedure. 


114112 

Ikfcib 


(2.156) 


Algorithm  B 

Comment  [Define  the  computational  grid.] 

Create  Af  =  2"*  suLjntervals  on  [a,  c]  by  choosing  a  sequence  of  boundary  points  61, 62*  •  •  • , 
with  6’  =  a  and  b^+i  =  c.  Choose  the  number  p  of  Chebyshev  nodes  on  each 
interval  flj"  =  [6i,  6«+i]  for  t  =  1, . . .,  M.  Determine  the  locations  of  the  scaled  Chebyshev 
nodes  ,r} , . . .  on  each  interval  B^' ,  and  use  the  initial  approximations  $oi  ^0  to 
evaluate  the  initial  approximations  $,<7  at  these  nodes,  obtaining  Choose 

tolerance  c. 


Step  1. 

Comment  [Use  Algorithm  A  to  compute  Newton  iterates  $4,  obtaining  the  solution  $  of  equation 
(2.5).] 

repeat 

(1)  Set  4  =  O’  =  <T. 

(2)  Evaluate  the  functions  $1,  g  at  each  of  the  scaled  Chebyshev  nodes  rl,  r?, . . .  ,7f  on 
each  interval  obtaining 

(3)  Apply  Algorithm  A  to  the  discretized  form  of  (2.50),  obtaining  6. 

(4)  Set  <T  =  d  +  6. 

(5)  Compute  the  solution  $  of  equation  (2.49)  from  the  values  of  <t,  using  composite 
Chebyshev  quadrature  (see  Step  6  of  Algorithm  A  and  Remark  2.11  above). 

until  |bll2/|ld||2  <  6 
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2.4  Numerical  Results 

FORTRAN  programs  have  been  written  implementing  the  algorithms  described  in  the  pre¬ 
ceding  section.  In  this  section,  we  discuss  several  details  of  our  implementation,  and  demon¬ 
strate  the  performance  of  the  scheme  with  numerical  examples. 

The  following  technical  details  of  our  implementation  appear  to  be  worth  mentioning. 

1.  The  algorithms  described  in  the  preceding  section  require  that  the  number  M  of  ele¬ 
mentary  subintervals  on  the  interval  [a,  c]  be  a  power  of  2.  Clearly,  this  is  not  an  essential 
limitation  and  it  can  be  removed  by  simple  bookkeeping  changes.  In  the  version  of  the 
algorithms  used  for  numerical  experiments,  these  changes  were  made. 

2.  Algorithm  A  depends  for  its  stability  on  the  equations  (2.140),  (2.141)  having  unique 
solutions  for  all  subintervals  B-  (/  =  0, 1, . . . ,  M,  i  =  1, .  • . ,  2^),  while  Algorithm  B  depends 
on  (2.152),  (2.153)  having  unique  solutions  for  all  subintervals  and  for  all  Newton  it¬ 
erates  k.  It  is  easy  to  construct  examples  for  which  these  conditions  are  violated,  even 
though  equation  (2.57)  or  equation  (2.32)  has  a  unique  solution.  In  such  cases,  a  different 
subdivision  of  the  interval  [a,  c]  can  be  attempted,  such  that  none  of  the  subintervals 

of  the  new  subdivision  coincides  with  an  interval  of  the  original  one.  This  procedure  can 
be  viewed  as  a  form  of  pivoting,  and  it  is  easy  to  show  that  it  is  always  possible  to  make 
it  work.  It  has  not  been  implemented  at  this  point,  and  we  have  not  so  far  encountered  a 
need  for  it. 

3.  We  have,  however,  implemented  a  crude  scheme  for  detecting  high  condition  numbers  in 
the  algorithms.  These  can  occur  in  two  places:  in  the  solution  of  the  linear  systems  on  each 
of  the  finest  level  subintervals  (Step  1  of  Algorithm  A),  and  while  computing  coefficients 
Ai ,  A2  defined  by  (2.90),  (2.91)  used  when  merging  solutions  on  two  consecutive  subintervals 
(Step  3  of  Algorithm  A).  In  both  cases,  the  condition  number  of  the  system  being  solved  is 
estimated  in  the  process  of  solution  (we  use  a  standard  UNPACK  routine),  and  the  largest 
of  these  is  returned  to  the  user.  When  an  extremely  large  condition  number  is  detected 
by  the  UNPACK  routine,  the  resulting  solution  of  the  original  ODE  should  be  viewed  as 
suspect.  It  is  easy  to  show  that  when  the  differential  operator  is  positive  definite,  this 
cannot  happen.  A  more  complete  treatment  of  this  subject  requires  further  study. 

4.  In  the  upward  sweep  (Step  3)  of  Algorithm  A,  we  evaluate  the  matrices  q*/,q*r  for 

all  intervals  Bij  and  use  these  matrices  to  evaluate  the  vectors  the  vectors  A'’^, 

and,  finally,  the  solution  a  of  the  integrzJ  equation  (2.57).  But  the  matrices  a'/,aR  do  not 
depend  on  the  right-hajid  side  /  of  equation  (2.57),  and  it  is  easy  to  see  that  their  evaluation 
accounts  for  more  than  90%  of  the  work.  Therefore,  whenever  the  equation  (2.57)  has  to  be 
solved  with  multiple  right-hand  sides,  we  can  precompute  the  matrices  a'£‘,a'R  and  store 
them,  saving  90%  of  the  cost  of  the  evaluation  of  subsequent  solutions. 

The  algorithms  of  this  chapter  have  been  applied  to  a  variety  of  problems.  Five  experi¬ 
ments  are  described  below,  and  their  results  are  summarized  in  Tables  2.1-2.13. 

Tables  2.1-2.11  are  associated  with  examples  for  which  analytic  solutions  are  available. 
In  each  of  these  tables,  the  first  column  contains  the  total  number  N  of  nodes  in  the 
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discretization  of  the  interval  [a,c].  The  second  column  contains  the  relative  L^  error  of 
the  numerical  solution  compared  with  the  analytically  obtained  one  at  5000  equispaced 
points  within  the  interval  [a,c],  where  Chebyshev  interpolation  has  been  used  to  evaluate 
the  numerical  solution  at  each  of  the  5000  points.  The  third  column  contains  the  maximum 
absolute  error  obtained  at  any  of  the  5000  points.  The  fourth  column  contains  the  CPU 
time  required  to  solve  the  problem,  excluding  the  time  used  to  evaluate  the  solution  at 
5000  equispaced  points,  where  in  all  cases  the  times  are  given  for  a  SUN  SPARCstation  1 
computer.  Tables  2.9-2.11,  associated  with  a  nonlinear  example,  have  in  addition  a  fifth 
column  which  contains  the  number  of  Newton  steps  taken  before  the  stopping  criterion 
(2.157)  has  been  satisfied,  with  e  =  10“^°. 

Tables  2.12-2.13  are  associated  with  an  example  for  which  we  did  not  have  analytic 
solutions.  In  this  example,  we  compare  each  numerical  solution  with  p  Chebyshev  nodes 
and  n  subintervals  against  the  solution  with  p  Chebyshev  nodes  and  2-n  subintervals.  In  each 
of  these  tables,  the  first  column  contains  the  total  number  N  of  nodes  in  the  discretization 
of  the  interval  [a,c].  The  second  column  contains  the  relative  L^  error  of  the  numerical 
solution  as  compared  with  the  numerical  solution  with  twice  the  number  of  subintervals, 
where  the  comparison  is  made  at  each  of  5000  equispaced  points  in  the  interval  [a,  c],  and 
where  Chebyshev  interpolation  has  been  used  to  evaluate  the  numerical  solution  at  each  of 
the  5000  points.  The  third  column  contains  the  maximum  absolute  error  obtained  at  any  of 
the  5000  points.  The  fourth  column  contains  the  SPARCstation  CPU  time  required  to  solve 
the  problem,  excluding  the  time  used  to  evaluate  the  solution  at  5000  equispace  points. 

Remark  2.12  In  Example  2.3  below,  we  solve  a  system  of  boundary  value  problems  of 
order  2;  and  in  Example  2.5,  we  solve  a  problem  of  order  4.  In  both  cases,  the  problems 
were  reduced  to  canonical  first  order  systems  (see,  for  example,  [13]),  with  the  latter  solved 
by  means  of  algorithms  A  or  B  of  the  preceding  section,  as  appropriate.  □ 

Example  2.1  This  example  is  taken  from  [6],  where  it  is  introduced  as  a  stiff  problem. 
The  equation  to  be  solved  is  given  by  the  formulae 


<^'i(a:)  -  998- <?l)i(i)  -  1998- <?!>2(x)  =  2  •  x,  (2.158) 

0^(x)  +  999  •<?^i(x)+  1999- <^(x)  =  x,  (2.159) 

subject  to  the  boundary  conditions 

<^i(0)  -  1,  (2.160) 

<;^2(1)  =  f  5-e-'®°”  +  .004-(.999  +  .001-e“^°°°).  (2.161) 

We  use  the  results  of  Theorem  2.3  to  reduce  the  first  order  system  (2.158)-(2.161)  to  one 

subject  to  homogeneous  boundary  conditioi 

0,(0)  =  0,  (2.162) 

02(1)  =  0.  (2.1G3) 
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Table  2.1:  Numerical  results  for  Example  2.1,  p  =  8. 


n 

E^(¥) 

E°°(*) 

t  (sec.) 

16 

0.962  xlO" 

0.685  xlO^ 

0.150  xlO" 

32 

0.244  xlO* 

0.216  xlO^ 

0.300  X 10° 

64 

0.700  xlO"* 

0.272  xlO^ 

0-560  xl0° 

128 

0.255  xlO-^ 

0.125  xlO^ 

0.108  xlO^ 

256 

0.667  xlO-2 

0.342  xl0° 

0.214  xlO^ 

512 

0.805  xlO-3 

0.536  xlO"* 

0.428  xlO^ 

1024 

0.330  xlO--* 

0.291  xlO-2 

0.846  xlO^ 

2048 

0.482  xlO-® 

0.529  xlO"”* 

0.168  XxG^ 

4096 

0.316  xlO-* 

0.476  xlO-® 

0.337  xlO^ 

8192 

0.112  xlO-^o 

0.170  xlO"* 

0.670  xlO^ 

16384 

0.115  xlO-^^ 

0.113  xl0-^° 

0.137  xlO^ 

We  apply  Algorithm  A  to  this  system  using  equispaced  subintervals,  with  the  number  of 
Chebyshev  nodes  p  =  8,16,24.  For  this  experiment,  the  background  Green’s  function  is 
chosen  to  correspond  to  the  equation 

$'(x)  =  0,  (2.164) 

subject  to  boundary  conditions  (2.162)-(2.163).  The  results  of  this  experiment  are  presented 
in  Tables  2. 1-2.3. 

Example  2.2  We  solve  the  problem  (2.158)-(2.161)  defined  in  Example  2.1,  but  using 
an  alternate  division  of  the  subintervals.  Since  the  solution  of  this  problem  has  a  fairly 
sharp  boundary  layer  near  the  left  end  of  the  interval  [0, 1],  we  construct  the  intervals 
via  the  formula 

6i  =  0, 

bi  =  f-j  for  i  =  2,. . ., M  +  1  ,  (2.165) 

so  that  they  become  progressively  smaller  near  the  left  end  of  the  interval  [0,1].  As  in 
Example  2.1,  we  reduce  the  problem  (2.158)-(2.161)  to  a  first  order  system  subject  to 
homogeneous  boundary  conditions  (2.162)-(2.163).  Algorithm  A  has  been  appbed  to  this 
problem  using  the  Green’s  function  corresponding  to  the  equation  (2.164)  subject  to  bound¬ 
ary  conditions  (2.162)-(2.163),  and  with  the  number  of  Chebyshev  nodes  p  =  16  and  24. 
The  results  of  this  experiment  appear  in  Tables  2.4-2.5,  and  are  most  satisfactory. 

Example  2.3  The  purpose  of  this  example  is  to  demonstrate  the  performance  of  the 
method  when  the  coefficient  p  of  the  equation  (2.2)  is  singular  at  the  ends  of  its  inter¬ 
val  of  definition,  while  the  particular  solution  being  sought  is  smooth.  We  solve  the  Bessel 
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Table  2.2:  Numerical  results  for  Example  2.1,  p  =  16. 


n 

F(¥) 

t  (sec.) 

16 

0.567  xlO" 

0.436  xlO^ 

0.310  xlO“ 

32 

0.251  xlO^ 

0.192  xlO^ 

0.540  xlO° 

64 

0.340  xlO-^ 

0.139  xlO^ 

0.950  xlO° 

128 

0.744  xlO-2 

0.344  xlO° 

0.179  xlO^ 

256 

0.798  xlO'^ 

0.389  xlO“^ 

0.345  xlO* 

512 

0.164  xlO-'* 

0.930  xlO-2 

0.686  xlO^ 

1024 

0.364  xlO"’’ 

0.243  xlO-5 

0.137  xlO^ 

2048 

0.992  xl0-“ 

0.817  xlO"^ 

0.274  xlO^ 

4096 

0.942  xlO"^^ 

0.776  xl0~^2 

0.532  xlO^ 

8192 

0.214  xlO-^2 

0.164  xlO-" 

0.107  xlO^ 

Table  2.3:  Numerical  results  for  Example  2.1,  p  =  24. 


n 

t  (sec.) 

24 

0.511  xlO" 

0.373  xlO^ 

0.690  xlO*^ 

48 

0.251  Xl0° 

0.244  xlO^ 

0.116  xlO^ 

96 

0.591  xlO-2 

0.240  xlO° 

0.209  xlO^ 

192 

0.496  xlO"^ 

0.214  xlO"^ 

0.387  xlO^ 

384 

0.546  xlO-5 

0.258  xlO-3 

0.765  xlO^ 

768 

0.274  xlO"® 

0.129  xlO"® 

0.147  xiO^ 

1536 

0.105  xlO-^2 

0.335  xlO-” 

0.295  xlO^ 

3072 

0.663  X  10-^3 

0.624  X  10-^2 

0.578  xlO^ 

Table  2.4:  Numerical  results  for  Example  2.2,  p  =  16. 


n 

W{¥) 

E°°($) 

i  (sec.) 

16 

0.567  xlO“ 

0.436  xlO* 

0.310  xlO^i 

32 

0.251  xlO^ 

0.192  Xl02 

0.540  xlO° 

64 

0.790  xlO-2 

0.360  X 10° 

0.950  xlO° 

128 

0.992  xlO-^^ 

0.818  xlO-° 

0.177  xlOi 

256 

0.244  X 10-12 

0.261  X 10-11 

0.352  xlOi 

512 

0.243  X 10-12 

0.261  X 10-11 

0.681  xlOi 
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Table  2.5:  Numerical  results  for  Example  2.2,  p  —  24. 


n 

£^(¥) 

£oo($) 

t  (sec.) 

24 

0.511  xl0“ 

0.373  xlQi 

0.710  xl0“ 

48 

0.251  xlO® 

0.244  xlQi 

0.117  xlOi 

96 

0.496  xl0~^ 

0.214  xlO-i 

0.209  xlOi 

192 

0.293  X 10-12 

0.227  X 10-11 

0.387  xlOi 

384 

0.294  X 10-12 

0.227  X 10-11 

0.758  xlOi 

equation  system 

^ - -  ■  4>v{x)  +  -  •  =  0,  (2.166) 

X 

2  2  1 

+  -  •  4>u-2{x)  =  0,  (2.167) 

X^  X 

-  -  •  4>u-i{x)  +  +  5  •  n  -  6  ^  ^2.168) 

X  X 

(see,  for  example,  [1])  on  the  interval  [0,600]  with  the  boundary  conditions 

<^.^(0)  =  «^^_i(0)  =  <^^_2(0)  =  0,  (2.169) 

<?i(,(600)  =  0.030598170290372796,  (2.170) 

0t-i(6OO)  =  0.015416721257492013,  (2.171) 

</»!.-2(600)  =  -0.025526503991812874,  (2.172) 


and  u  =  100.  The  difficulty  of  this  problem  is  due  to  the  fact  that  the  two  linearly  inde¬ 
pendent  solutions  to  each  of  equations  (2.166),  (2.167),  (2.168)  are  J„(r),  Tl/(x); 
yi,_i(r);  and  /i,_2(i),  yi,_2(x),  respectively,  (Bessel  functions  of  the  first  and  the  second 
kinds).  As  is  well  known,  7^(i),  and  J^-2  behave  in  the  vicinity  of  zero  like  x'',x‘'~^, 
and  x^~^,  respectively,  while  yj,(a;),yj,_i,  and  yi,_2  behave  like  x~‘' ,x~^‘'~^\  and  x~^‘'~^\ 
respectively;  most  methods  have  trouble  finding  the  decaying  solution.  In  addition,  this  is 
a  fairly  large-scale  calculation,  since  the  the  solution  to  (2.166)-(2.168)  contains  almost  100 
wavelengths  in  the  interval  [0,600]. 

We  reduce  the  problem  (2.166)-(2,168)  to  a  first  order  system,  and  apply  Algorithm  A  to 
this  system  using  equispaced  subintervals,  with  the  number  of  Chebyshev  nodes  p  =  16,20 
and  24.  For  this  experiment,  the  background  Green’s  function  is  chosen  to  correspond  to 
the  equation 

$'(x)  =  0,  (2.173) 

subject  to  boundary  conditions  (2.169)-(2.172).  The  results  of  this  experiment  are  presented 
in  Tables  2. 6-2. 8. 
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Table  2.6:  Numerical  results  for  Example  2.3,  p  =  16. 


n 

E\^) 

t  (sec.) 

128 

0.133  xlO^ 

0.415  xlO'^ 

0.193  xlO^ 

256 

0.958  xl0° 

0.110  xl0° 

0.384  xlO^ 

512 

0.758  xlO-2 

0.779  xlO-3 

0.764  xlO^ 

1024 

0.203  xlO-5 

0.213  xlO-® 

0.152  xlO^ 

2048 

0.632  xl0-i° 

0.661  xlO-" 

0.301  xlO^ 

4096 

0.823  xlO-^^ 

0.126  xlO"” 

0.609  xlO^ 

Table  2.7:  Numerical  results  for  Example  2.3,  p  =  20. 


n 

F(¥) 

£;oo($) 

t  (sec.) 

160 

0.185  xlO^ 

0.558  xlO“ 

0.330  xlO'-^ 

320 

0.216  xlO' 

0.392  xlO° 

0.655  xl02 

640 

0.373  xlO-^ 

0.388  xlO"® 

0.130  xlO^ 

1280 

0.493  xlO"^ 

0.522  xlO-^° 

0.260  xlO^ 

2560 

0.111  xlO-^2 

0.220  xlO'^^ 

0.520  xlO^ 

5120 

0.142  xlO-'° 

0.186  xlO-" 

0.104  xlO"* 

Table  2.8:  Numerical  results  for  Example  2.3,  p  =  24. 


n 

E\^) 

t  (sec.) 

96 

0.125  xlO^ 

0.431  xlO*^ 

0.262  xlO^ 

192 

0.133  xlO^ 

0.277  xlO° 

0.527  xlO^ 

384 

0.718  xlO-i 

0.797  xlO-2 

0.104  xlO^ 

768 

0.728  xlO"'^ 

0.765  xlO"® 

0.206  xlO^ 

1536 

0.308  X  10-^2 

0.555  xlO-^^ 

0.409  xlO^ 

3072 

0.551  xl0-'2 

0.397  xlO-^3 

0.818  xlO^ 
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Example  2.4  Wa  consider  a  system  of  Jacobian  elliptic  functions  sn,  cn,  dn  :  [0,  lO  A']  — >  R 
(see,  for  example,  [1])  which  are  solutions  to  the  equations 


sn'{x)  =  cn(ar)  •  dn(i), 
cn'{x)  =  -sn(x)  •  dn{x), 
dn'(x)  =  —Tn-sn{x)-cn{x), 


with  m  =  ^  in  our  experiments,  subject  to  the  boundary  conditions 


(2.174) 

(2.175) 

(2.176) 


sn(0)  =  0, 
cn(0)  =  1, 
dn(40-K)  =  1, 


with  K  given  by  the  expression 

, 

Jo  vl  —  m  •  sin^  9 


(2.177) 

(2.178) 

(2.179) 


(2.180) 


VVe  use  for  an  initial  guess  the  solution  to  (2.176),  (2.179)  for  m  =  0,  which  is  defined 
by  the  formulae 


sn(a:)  =  sin  ^  ,  (2.181) 

cn(i)  =  cos^2^’*)’  (2.182) 

dn{x)  =  1.  (2.183) 


We  use  the  results  of  Theorem  2.4  to  reduce  thp  system  (2.174)-(2.176)  to  one  subject 
to  the  homogeneous  boundary  conditions 


sn(0)  =  0, 
cn(0)  =  0, 
<in(40ii')  =  0, 


(2.184) 

(2.185) 

(2.186) 


and  then  apply  Algorithm  B  to  this  system  using  equispaced  subintervals,  with  the  number 
of  Chebyshev  nodes  p  =  8, 16  and  32.  For  this  experiment,  the  background  Green’s  function 
is  chosen  to  correspond  to  the  equation 


$'(i)  =  0, 

subject  to  boundary  conditions  (2.184)-(2.186).  The  results  of  this  experiment  are  presented 
in  Tables  2.9-2.11. 

Example  2.5  This  example  is  taken  from  [28].  Its  purpose  is  to  demonstrate  the  perfor¬ 
mance  of  the  method  when  the  equation  to  be  solved  contains  fourth  order  derivatives.  The 


38 


CHAPTER  2.  TWO-POINT  BOUNDa^i  .  v^ALUE  PROBLEMS 


Table  2.9:  Numerical  results  for  Example  2.4,  p  =  8. 


n 

W{¥) 

t  (sec.) 

Steps 

128 

0.211  xlO-i 

0.551  xlO-* 

0.167  xlO^ 

9 

256 

0.158  xlO-2 

0.416  xlO-2 

0.258  xlO^ 

7 

512 

0.849  xlO-® 

0.226  xl0-‘‘ 

0.439  xlO^ 

6 

1024 

0.106  xlO-^ 

0.330  xlO-^ 

0.883  xlO^ 

6 

2048 

0.313  xl0-^° 

0.984  xl0-*“ 

0.175  xlO^ 

6 

4096 

0.147  xlO-^2 

0.469  X  10-^2 

0.348  xlO^ 

6 

Table  2.10:  Numerical  results  for  Example  2.4,  p  =  16. 


n 

Wi¥) 

t  (sec.) 

Steps 

64 

0.162  xl0° 

0.505  xlO° 

0.177  xlO^ 

10 

128 

0.422  xlO-‘ 

0.108  xlO® 

0.246  xlO^ 

7 

256 

0.181  xlO-^ 

0.476  xlO'^ 

0.409  xlO^ 

6 

512 

0.441  xlO-^ 

0.120  xlO"® 

0.851  xlO^ 

6 

1024 

0.425  X  10-^2 

0.125  xlO"” 

0.164  xlO^ 

6 

2048 

0.115  xlO-12 

0.317  xlO-^2 

0.324  xlO^ 

6 

4096 

0.569  xlO"^^ 

0.152  xlO"^^ 

0.653  xlO^ 

6 

Table  2.11:  Numerical  results  for  Example  2.4,  p  =  32. 


n 

£^(¥) 

t  (sec.) 

Steps 

256 

0.148  xlO-^ 

0.388  xlO'^ 

0.163  xlO^ 

6 

512 

0.124  xlO-^ 

0.395  xlO-^ 

0.322  xlO^ 

6 

1024 

0.450  xl0-'2 

0.119  xlO-" 

0.632  xlO^ 

6 

2048 

0.266  x  10-^2 

0.703  xl0->2 

0.126  xlO" 

6 
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deflection  of  a  beam  under  a  uniform  load  q,  with  the  beam  built  in 
and  simply  supported  at  the  right  end,  is  given  by  the  formula 

at  the  left  end  (x  =  0) 

y""(x)  +  ;g^-sW  =  :E^. 

(2.1S7) 

subject  to  the  boundary  conditions 

o 

II 

II 

II 

o 

"»» 

1! 

o 

(2.188) 

with  k  the  force  per  unit  deflection  per  unit  length  of  beam,  and  E 
of  the  beam.  The  values  of  the  constants  used  are 

•  /  the  flexural  rigidity 

L  =  1.2  X  10^  in.. 

(2.189) 

k  —  2.604  X  10^  psi. 

(2.190) 

q  =  4.34  X  10'’  Ibs/in., 

(2.191) 

E  =  3.0  X  lO’^  psi. 

(2.192) 

/  =  3.0  X  10^  in.^ 

(2.193) 

(see  [28],  p.  174).  The  L^  norm  of  y,  the  solution  to  (2.188),  is  approximately  10®  times 
larger  than  the  L^  norm  of  y"".  Combined  with  the  high  number  of  derivatives  in  (2.187), 
this  tends  to  present  difficulties  for  finite  difference  methods.  We  reduce  the  problem  (2.187), 
(2.188)  to  a  first  order  system,  and  then  use  the  results  of  Theorem  2.9  to  express  $  by  the 
formula 

$(x)  =  'j»(x)  •  r(i), 

with  ^(x)  :  [0,  L]  -*  L(R‘‘^‘‘)  given  by  the  formula 


/  L~x 

L 

0 

X 

L 

0  \ 

0 

1 

0 

0 

0 

0 

0 

1 

V  -1 

0 

1 

and  r  ;  [a,c]  — ►  R"  the  solution  to  the  equation 

r(x)  +  »lr-»(x)  ■  (^r'(x)  +  p(x)  •  ^'(x))  •  r(x)  =  ^-^(x)  •  /(x). 


subject  to  boundary  conditions 


/ 

1 

0 

0 

0  \ 

{  ® 

0 

0 

0  \ 

0 

1 

0 

0 

■  9(a)  •r(a)  + 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

\ 

0 

0 

0 

0 ) 

1 1 

0 

0 

0/ 

■  «'(c)-r(c)  =  o. 


with  p  :  [0.  L]  -*  L(R'*’^'*)  defined  by  the  formula 


/  0 

-1 

0 

0 

\ 

0 

0 

-1 

0 

0 

0 

0 

-1 

1  JS— 

\  El 

0 

0 

0 

) 

(2.194) 


(2.195) 


(2.196) 


p(x)  = 


(2.197) 
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Table  2.12:  Numerical  results  for  Example  2.5,  p  =  4. 


n 

£^(¥) 

£“($) 

t  (sec.) 

4 

0.115  xlO" 

0.672  xlO"^ 

0.110  xlO** 

8 

0.755  xlO-2 

0.506  xlO-2 

0.220  xl0° 

16 

0.474  xlO-^ 

0.369  xlO-3 

0.400  xlO^ 

32 

0.269  xlO"'* 

0.249  xlO-** 

0.780  xl0° 

64 

0.185  xlO-5 

0.162  xlO-5 

0.153  xlO^ 

128 

0.117  xlO-® 

0.103  xlO-® 

0.303  xlO^ 

256 

0.891  xlO-* 

0.652  xlO-® 

0.603  xlO^ 

512 

0.306  xlO"* 

0.170  xlO"* 

0.120  xlO^ 

Table  2.13:  Numerical  results  for  Example  2.5,  p  =  8. 


n 

P(¥) 

E°°i4>) 

t  (sec.) 

8 

16 

32 

64 

0.130  xlO-^ 
0.792  xlO-® 
0.465  xlO-* 
0.916  xl0-« 

0.816  xlO-*^ 
0.510  xlO-® 
0.235  xlO'* 
0.463  xlO-® 

0.280  xl0“ 
0.530  xl0° 
0.105  xlO^ 
0.201  xlO^ 

and  /  :  [0,  L]  -*  R"*  defined  via  the  formula 


/  0  \ 
0 
0 

^  'ii  > 


(2.198) 


We  apply  Algorithm  A  to  this  problem  using  equispace  subintervals,  with  the  number 
of  Chebyshev  nodes  p  =  4  and  8.  For  this  experiment  the  background  Green’s  function  is 
chosen  to  correspond  to  the  equation 


$'(a:)  =  0, 

subject  to  boundary  conditions  (2.188).  The  results  of  this  experiment  are  presented  in 
Tables  2.12-2.13. 

The  following  observations  can  be  made  from  Tables  2.1-2.13,  and  are  corroborated  by 
our  more  extensive  experiments. 

1.  The  practical  convergence  rate  of  the  method  is  consistent  with  the  theoretical  one.  For 
larger  p.  the  exact  numerical  verification  of  the  order  of  convergence  tends  to  be  difficult, 
since  the  precision  of  calculations  is  exhausted  before  the  behavior  of  the  scheme  becomes 
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asymptotic.  However,  this  is  often  encountered  when  dealing  with  rapidly  convergent  algo¬ 
rithms. 

2.  For  small-scale  problems  (such  as  in  Example  2.5)  and  large  p,  the  algorithm  produces 
essentially  exact  results  with  a  small  number  of  nodes.  For  large-scale  problems,  double 
precision  accuracy  is  achieved  at  approximately  20  nodes  per  wavelength  with  p  =  20,  at 
12  nodes  per  wavelength  with  p  =  24,  and  at  10  nodes  per  wavelength  with  p  =  32.  The 
optimal  timings  are  achieved  at  p  between  24  and  32  (provided  that  about  10-12  digits  of 
accuracy  are  desired). 

3.  The  condition  number  of  a  Nystrom  discretization  of  a  second  kind  integral  equation 
is  asymptotically  bounded,  and  our  results  reflect  this  fact.  The  relatively  poor  accuracy 
(8-11  digits)  obtained  in  Example  2.5  is  due  to  the  ill-conditioning  of  the  original  ODE,  as 
opposed  to  that  of  the  numerical  scheme  used. 

4.  The  algorithm  is  completely  indifferent  to  the  stiffness  near  the  left  end  of  the  interval 
[0, 1]  of  equations  (2.158),  (2.159)  in  Examples  2.1-2.2. 

5.  It  is  easy  to  use  the  algorithm  in  an  adaptive  manner,  as  demonstrated  in  Example 
2.2.  However,  a  fully  adaptive  version  of  the  scheme  has  not  been  implemented.  The 
intervals  in  Example  2.2  were  provided  by  the  calling  program  (as  opposed  to  having 
been  constructed  by  the  algorithm  itself). 

6.  If  the  function  p  ;  [a,  c]  L(R'**'”)  given  by  (2.1)  is  singular  in  the  interval  [a,c],  then 

the  choice  of  a  background  Green’s  function  can  dramatically  affect  the  numerical  results 
(see  Remark  2.8).  When  p  is  not  singular  in  the  interval  [a,c],  the  numerical  advantages 
of  one  background  Green’s  function  over  another  are  usually  minor.  However,  using  the 
Green’s  function  given  in  Lemma  2.2  results  in  a  slightly  faster  algorithm.  This  is  because 
this  Green’s  function  is  constant  in  each  of  the  intervals  (i  <  t),  {x  >  t),  which  provides 
in  Step  2  of  Algorithm  A  faster  evaluation  of  the  matrices  given  by  equations 

(2.143)-(2.144)  and  vectors  given  by  equations  (2.145)-(2.146),  and  provides  in 

Step  6  a  faster  evaluation  of  the  solution  $  of  equation  (2.27). 

7.  The  algorithm  caji  solve  systems  of  high  order  equations  with  no  numerical  difficulty,  as 
demonstrated  by  Example  2.5. 
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Chapter  3 

One-Dimensional  Integral 
Equations 


3.1  Introduction 

Id  tliis  chapter,  we  consider  the  problem  of  determining  a  function  a  :  [a,c]  —>■  R  which 
satisfies  the  integral  equation 

A  •  <T(a:)  +  /  k{x,t)  ■  <r{t)dt  =  f(x),  (3.1) 

Ja 

where  the  free  term  /  :  [a,  c]  — >  R  and  the  kernel  k  :  [o,  c]  x  [o,  c]  -+  R  are  known  functions, 
and  0  <  A  <  1  (so  that  (3.1)  is  a  first  kind  equation  when  A  =  0,  and  is  a  second  kind 
equation  when  A  =  1).  We  assume  that  the  kernel  k  contains  a  singularity  of  the  form 


s(x-t)  =  logii-tl. 

(3.2) 

or  of  the  form 

s(i  - 1)  =  ji  -  tr,  (0  <  iqi  <  1), 

(3.3) 

or  of  the  form 

s(x  —  t)  =  — — . 

^  ’  x-t 

(3.4) 

When  s  is  given  by  (3.2)  or  (3.3),  equation  (3.1)  is  a  weakly  singular  integral  equation;  for 
a  singularity  given  by  (3.4),  equation  (3.1)  is  a  Cauchy  integral  equation,  and  its  integral  is 
evaluated  in  the  principal  value  sense. 

For  purposes  of  analysis,  equations  of  the  form  (3.1)  are  usually  divided  into  three 
classes: 

(1)  Second  kind  integral  equations  with  singularities  of  the  form  (3.2)  or  (3.3). 

(2)  First  and  second  kind  integral  equations  with  singularities  of  the  form  (3.4). 

(3)  First  kind  integral  equations  with  second  singularities  of  the  form  (3.2)  or  (3.3). 
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Generally  speaking,  integral  operators  for  any  of  these  classes  have  eigenfunctions  with 
end-point  singularities,  so  that  an  integral  equation  will  have  a  smooth  solution  only  if  the 
solution  is  orthogonal  to  such  eigenfunctions.  Each  class  of  equations  also  has  its  own  set  of 
theorems  governing  the  existence  and  uniqueness  of  solutions.  When  second  kind  integral 
equations  with  weak  singularities  are  considered,  Fredholm’s  theorems  govern  the  existence 
and  uniqueness  of  solutions,  and  the  linear  systems  arising  from  discretization  axe  generally 
well-conditioned.  When  first  and  second  kind  integral  equations  with  Cauchy  singularities 
are  considered,  Fredholm’s  theorems  apply  in  an  extended  sense  (see  Remark  3.3),  and 
the  discretized  linear  systems  are  also  generally  well-conditioned.  Fredholm’s  theorems  do 
not  apply  to  first  kind  integral  equations  with  weak  singularities,  although  some  first  kind 
integral  equations  are  known  to  have  unique  solutions.  First  kind  integral  equations  also 
yield  ill-conditioned  linear  systems:  the  condition  number  of  the  discretized  system  is  at 
lecLSt  0(N)  where  N  is  the  dimension  of  the  problem. 

When  an  integral  equation  is  discretized  using  a  Nystrom  scheme,  the  order  of  conver¬ 
gence  of  the  method  is  equivalent  to  the  order  of  convergence  of  the  underlying  quadrature 
formula.  Standard  quadrature  formulas  yield  extremely  poor  convergence  for  the  integral 
equations  considered  in  this  chapter,  due  to  the  kernel  singularities  in  these  equations.  Spe¬ 
cial  quadrature  formulae  have  been  developed  for  these  types  of  integral  equations.  When 
the  kernel  has  a  Cauchy  singularity,  the  methods  presented  in  [8]  and  [9]  yield  quadrature 
formulas  with  superalgebraic  convergence.  When  the  kernel  contains  a  weak  singularity, 
the  quadrature  formulas  developed  in  [27],  [3]  yield  up  to  eighth-order  convergence. 

In  this  chapter,  we  construct  two  rapidly  convergent,  order  0(N)  algorithms  for  the 
direct  solution  of  first  and  second  kind  integral  equations  containing  weak  singularities  or 
Cauchy  singularities.  The  first  algorithm  is  designed  for  integral  equations  with  end-point 
singularities  in  the  solution,  while  the  second  algorithm  is  designed  for  integral  equations 
with  smooth  solutions.  We  extend  the  observations  of  [5]  to  construct  sparse  representations 
of  the  integral  operators;  we  then  extend  the  techniques  of  [22]  and  [29]  to  construct  fast, 
direct,  rapidly  convergent  solvers.  In  addition,  we  extend  the  methods  of  [27]  and  [3]  to  yield 
superaJgebraically  convergent  quadrature  formulas  for  equations  with  weak  singulcirities. 

The  plan  of  this  chapter  is  as  follows:  in  Section  3.2  we  review  the  relevant  properties 
and  tools  of  Chebyshev  approximation,  in  Section  3.3  we  apply  Chebyshev  analysis  to 
the  integral  operators  of  interest,  in  Section  3.4  we  develop  superaJgebraically  convergent 
quadrature  formula  for  integral  operators  with  weak  singularities;  in  Section  3.5  we  develop 
the  analytical  apparatus  for  the  algorithm  for  solutions  with  end-point  singularities;  in 
Section  3.6  we  develop  the  analytical  apparatus  for  the  algorithm  for  smooth  solutions; 
Sections  3.7  and  3.8  describes  the  numerical  schemes  for  the  algorithm  for  solutions  with  end¬ 
point  singularities  and  the  algorithm  for  smooth  solutions,  respectively;  finally,  in  Section 
3.9  we  illustrate  the  performance  of  the  algorithms  with  several  practical  examples. 
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3.2  Chebyshev  Approximation 

In  this  section,  we  define  operators  and  summarize  results  from  Chebyshev  approximation 
theory.  Most  of  the  results  are  classical  and  can  be  found,  for  example,  in  [16],  [18],  [20]. 
Much  of  our  discussion  follows  the  presentation  given  in  [22]. 


3.2.1  One-Dimensional  Chebyshev  Approximation 

Given  a  function  /  6  C*[-l,l],  we  define  the  vector  V  G  of  Chebyshev  coefficients  by 
the  expression 

y  p  (3.5) 

J-l  y/r^  ’  ^  ^ 

where  Ti  denotes  the  I'th  Chebyshev  polynomial.  The  function  /  can  therefore  be  represented 
by  the  expression 

(3.6) 

with  ►  L^[— 1, 1]  the  Chebyshev  interpolation  operator  given  by 

^ix)-V  =  f^Ti{x).Vi.  (3.7) 

f=0 


Let  t;  €  R'’  be  defined  by 

for  t  =  0, 1, . .  .,p  -  1,  and  let  rp  :  R'’ 


Ui  =  V,-, 

X^[-l,  1]  be  the  interpolation  operator  given  by 


p-i 

V’(x)- V  =  ^r,(z)-  ui. 

«=o 

We  denote  the  transpose  interpolation  operator  :  L^[— 1, 1] 

(^^(x)./).  =  v;-. 


(3.8) 

R'’  by  the  expression 

(3.9) 


for  i  =  1,. .  .,p,  where  VJ  is  given  by  (3.5). 

The  following  lemma  proves  that  a  Chebyshev  expansion  converges  rapidly  for  suffi¬ 
ciently  smooth  functions,  and  is  proved,  for  example,  in  [20]. 

Lemma  3.1  If  f  G  C*'[-l,  1],  u  G  R'’  is  given  by  (3.5)  for  i  =  0, 1, . .  .,p-  1,  and  rp  :W 
L^[-IA]  is  defined  by  (3.8),  then  for  any  x  G  [—1, 1], 

||/(x)-V’(x)-u||  =  0(^).  (3.10) 

Lemma  3.2  proves  that,  given  a  vector  /  G  R’’  representing  the  function  /  discretized 
at  the  roots  the  pth  Chebyshev  polynomial,  the  vector  v  G  R’’  given  by  (3.5)  for  i  = 
0,  l,...,p  -  1  may  be  obtained  from  /  via  a  discrete  cosine  transform.  A  proof  of  this 
lemma  may  be  found  in  [20]. 
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Lemma  3.2  Suppose  that  f  £  C^[0, 1],  v  €  R*’  is  given  by  (3.5)  for  i  =  0,1, . .  .,p  —  I,  and 
f  €  is  given  by  the  expression 

h  =  fiu),  (3.11) 

with  ti  the  ith  root  of  the  Chebyshev  polynomial  Tp{x).  Suppose  further  that  a  G  R^  is  given 
by  the  formula 

a=Cif),  (3.12) 

where  C  denotes  the  discrete  cosine  transform  of  dimension  p.  Then 

=  (3.13) 

Remark  3.1  Since  the  discrete  cosine  transform  of  a  function  may  be  obtained  via  the 
Fast  Fourier  Transform,  the  vector  of  Chebyshev  coefficients  a  may  be  obtained  from  the 

function  values  /  using  0(p  -  log  p)  arithmetic  operations.  □ 


While,  strictly  speaking,  Chebyshev  interpolation  is  defined  only  for  the  interval  [—1,1], 
we  can  define  an  interpolation  operator  t/’[a,c]  •  R’’  for  the  interval  [a,  c]  via  the 

formula 

^la,c]{x)  =  —  2  •  ^(0  +  ,  (3.14) 

so  that 

X  €  [a,c]  t  €  [-1, 1]  (3.15) 

(see,  for  example,  [15]).  Similarly,  we  define  the  transpose  operator  '•  L'^W,c]  — >  R’’ 
via  the  formula 

=  +  (3.16) 

Given  a  function  R?  we  denote  the  Chebyshev  coefficients  for  that  function 

by  that 

V’la, c]  • ’^[a.c]  (3.17) 

is  an  approximation  to  f[a,c]-  Given  the  p  roots  f,  of  the  Chebyshev  polynomial  Tp{t),  we 
define  roots  for  the  interval  [a,c]  via  the  formula 


t, 


(o.c] 


(c-q)  ,  ,  (q  +  c) 
2  '  2 


(3.18) 


for  i  =  0, 1, . .  .,p  -  1. 

Let  b  —  (a  +  c)/2  denote  the  midpoint  of  [a,c],  and  let  ^[a,b]  ■  R^  L^\a,b],  fh[b,c]  ■  R^ 
L^{b,c]  denote  the  interpolation  operators  for  the  intervals  [a,  6]  and  [6,  c],  respectively.  We 
consider  the  problem  of  obtaining  the  expansions  V[a,fc]  ^'^d  given  the  expansion  (.]. 
If  the  function  /  associated  with  is  sufficiently  smooth,  then  we  may  use  c]  to 
obtain  the  approximations  /(<![„  ^j)  and  /(fij^^j)-  Then,  (,]  is  obtained  from  /(<:[^(,])  via 
the  cosine  transform  described  in  Lemma  3.2,  and  similarly  is  obtained  from  /(tijj,,,])- 
Thus,  each  of  the  mappings  uja  c]  — ‘  ^[a,6]  3-nd  i’[3,c]  ^^(6.0]  is  obtained  via  the  application  of 
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an  interpolation  operator  and  a  cosine  transform  operator.  Let  Ca,Cb  '■  ^  R^  denote 

the  operators  which  perform  the  mappings  Wfa  ,.]  -+  and  ujo  ,.]  ^[6,c]?  respectively,  so 

that 


Ca  ■  ^[a,c]  —  ^[a,i>]»  (3.19) 

C'b  •  W[a,c]  =  V[(,,c]-  (3.20) 

Ca,Cb  are  nonsingular  (see,  for  example,  [5]);  we  now  define  the  inverse  operators  Cd,  Cu  : 
RP  RP  via  the  expressions 

Cd  =  C^\  (3.21) 

Cu  =  Cb^.  (3.22) 


3.2.2  Two-dimensional  Chebyshev  Approximation 

Given  a  function  K(x,t)  :  [—1,1]  X  [—1,1]  — R,  the  following  three  lemmas  provide  the 
methods  by  which  K  may  be  approximated  via  Chebyshev  approximation,  and  define  the 
circumstances  under  which  these  approximations  are  rapidly  convergent.  The  lemmas  are 
direct  consequences  of  Lemma  3.1. 


Lemma  3.3  Given  a  function  K{x,t)  :  [-1,1]  X  [-1,1]  -*•  R,  suppose  that  for  fixed  t, 
K{x,t)  €  C*'[-l,l].  Let  Ml  :  I^[-l,l]  x  I^[-l,l]  x  I^[-l,l]  be  defined  for  all 
t  €  [- 1 , 1]  via  the  formula 


Mi{K,t)i  =  1'^ 


Kix,t)-Ti{x) 

v^l  _  3,2 


dx, 


(3.23) 


so  that  for  each  t,  Mi{K,t)  yields  a  Chebyshev  expansion  in  x  which  approximates  K{x,t). 
Let  mi  :  L^[-l,l]  x  Z/^[— 1,1]  ^  R^  x  L^[—l,l]  be  given  by  the  expression 


mi{K,t)i  =  Mi{K,t)i,  (3.24) 

for  i  =  0,  l,...,p  —  1  and  for  all  t  £  [—1,1].  Finally,  let  ij}  :  R^  Z^[— 1,1]  be  the 
interpolation  operator  defined  by  (3.8).  Then,  for  any  x,t  £  [-1, 1], 

\\Kix,t)  -  ^(x)  •  mi{K,t)\\  =  O  .  (3.25) 


Lemma  3.4  Let  K{x,t)  :  [-1,1]  X  [-1,1]  — ►  R,  and  suppose  that  for  fixed  x,  K{x,t)  £ 
C^[-l,  1].  Let  M2  :  L^[-l,  1]  x  L^[-l,  l]  -♦  L^[-l,  1]  x  he  defined  for  all  x  £  [-1, 1]  via 
the  formula 


M, 


fK,x)i  =  1'^ 


Kix,t)-Ti{t) 


dt, 


(3.26) 


so  that  for  each  x,  M2{K,x)  yields  a  Chebyshev  expansion  in  t  which  approximates  K{x,t). 
Let  m2  :  L^[-l,  1]  x  L^[-l,  1]  — *  Z/^[— 1, 1]  x  R^  be  given  by  the  expression 


m2(A',x),  =  M2{K,x)i, 


(3.27) 
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for  i  =  0,l,...,p  —  1  and  for  all  x  6  [—1,1].  Finally,  let  :  £^[—1,1]  -+  R’’  be  the 
interpolation  operator  defined  by  (3.9).  Then,  for  any  x,t  £  [—1, 1], 

\\K{x,t)  -  m2(K,x) .  =  O  .  (3.28) 


Lemma  3.5  Given  a  function  K(x,t)  :  [—1,1]  X  [—1,1]  R,  suppose  that  for  fixed  x, 
K(x,t)  €  C*[-l,l].  Suppose  further  that  for  each  fixed  t,  K(x,t)  £  C^[— 1,1].  Let  M3  : 
X^[-l,  1]  X  X^[— 1, 1]  — >  X  P  denote  the  two-dimensional  Chebyshev  expansion  defined  by 
the  formula 


y/l^  ) 


Ti{x) 


dx. 


(3.29) 


for  each  i  y.  j  £  {0, 1, . . .}  x  {0, 1, .. .}.  Let  m3  :  £^[—1, 1]  X  £^[-1, 1]  — >  R^  X  R^  6e  given 
by  the  expression 

(3.30) 

for  all  i  =  0, 1, ...  ,p  -  1  and  for  all  j  =  0, 1, . .  .,p  —  1.  Finally,  let  V’  :  R’’  — >  ■^'^[“1?  1]> 
:  £^[-1,1]  -+  RP  6e  the  interpolation  operators  defined  by  (3.8),  (3.9),  respectively. 
Then,  for  any  x,t  £  [-1,1], 


\\K(x,  t)  -  V^(x) .  m^(K)  ■  ^^(t)||  =  0  . 


(3.31) 


Remark  3.2  Let  K{xi,ti)  ;  R”  — R"  €  L(R’‘’'")  denote  an  n  x  n  discretization  of  K, 
and  suppose  that  n  >  p.  Suppose  further  that  V’n  ^  R’’  -*  R”  is  the  operator  mapping  a 
Chebyshev  expansion  in  x  to  function  values  at  n  points  i,,  and  suppose  :  R”  — >  RP  is 
the  operator  mapping  function  values  at  n  points  in  U  to  Chebyshev  expansions  in  t.  Then, 
using  Lemmas  3.3-3.5,  we  may  represent  K  in  one  of  four  ways: 

Case  1:  If  for  fixed  i,  K{x,t)  £  C^[-l,l],  and  also  for  fixed  t,  K{x,t)  £  C^[— 1,1],  then 
by  Lemma  3.5  we  may  represent  K  by  the  formula 


K  =  rP.-m3{K)-rp'^,  (3.32) 

with  m3  defined  by  (3.30).  Thus,  K  can  be  approximated  by  m3,  and  merely 

used  to  obtain  values  of  K  at  specific  points  x,t.  We  have  reduced  the  representation  K 
to  the  p^  representation  m3. 

Case  2;  If  for  fixed  t,  K(x.t)  £  C*'[— 1, 1],  then  by  Lemma  3.3  we  may  represent  K  by  the 
formula 

K  =  ',P„-rhi{K,t),  (3.33) 

where  mi  :  Rp  — ‘  R"  denotes  the  operator  mi  given  by  (3.24),  discretized  at  the  n  points 
t,.  K  can  therefore  be  represented  by  mi  using  n  p  coefficients;  this  is  less  than  the 
coefficients  used  in  K ,  but  more  than  the  p^  coefficients  used  to  represent  m3  in  Case  1. 
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Case  3:  If  for  fixed  i,  K{x,t)  €  C^[— 1, 1],  then  by  Lemma  3.4  we  may  represent  K  by  the 
formiila 

k{x,t)  =  m2{K,x)-xl}'^,  (3.34) 

where  m2  :  R"  — ♦  denotes  the  operator  m2  given  by  (3.27),  discretized  at  the  n  points 

Xi.  The  number  of  coeIRcients  required  to  represent  K  using  m2  is  p  ■  n —  the  same  as  the 
number  required  in  Case  2,  but  more  than  the  number  required  in  Case  1. 

Case  4;  If  for  fixed  x,  K{x,t)  ^  C*[-l,l],  and  for  fixed  t,  K{x,t)  ^  C*'[— 1,1],  then 
we  cannot  use  Chebyshev  approximation  to  reduce  the  number  of  coefficients  needed  to 
represent  K.  We  require  coefficients  to  represent  K,  more  than  the  p  •  n  coefficients 
required  in  Cases  2  and  3,  and  more  than  the  p^  coefficients  required  in  Case  1.  □ 

3.3  Chebyshev  Approximation  for  Singular  Integral  Oper¬ 
ators 

In  this  section,  we  present  techniques  for  the  efficient  representation  of  integral  operators 
with  singularities,  with  the  result  that  a  dimension  N  discretization  of  an  integral  operator 
with  a  singularity  may  be  accurately  represented  using  only  0(N)  elements. 

We  assume  the  kernel  k  :  [a,  c]  X  [a,  c]  — ►  R  is  of  the  form 

A:(x,  t)  =  A:i(z,  t)  •  s(x  -t)  +  k2{x,  t),  (3.35) 

where  the  singularity  s  defined  by  one  of  (3.2)-(3.4),  and  ki,k2  :  [a,  c]  x  [o,  c]  R  6 
C^[a,c]  X  [a,c].  For  convenience,  we  define  the  operator  P  ;  X^[a,c]  — »  L^[a,c]  by 

P(ct)(i)  =  A  •  cr(x)  +  f  k{Xft)  ■  a{t)dt,  (3.36) 

Ja 

so  that  the  equation  (3.1)  assumes  the  form 


Pa  =  f.  (3.37) 

Remark  3.3  While  P  is  defined  as  an  operator  yielding  solutions  in  L^[a,  c],  the  existence 
and  uniqueness  of  solutions  in  X^[a,  c]  depends  on  the  class  of  the  integral  operator  under 
consideration.  If  P  is  a  second  kind  integral  operator  with  a  weak  singularity,  then  by  the 
Fredholm  Alternative,  either  there  exist  unique  solutions  a  €  X^[a,c],  or  the  homogeneous 
equation 

Pa  =  0  (3.38) 

has  a  nontrivial  solution. 

For  Cauchy  integral  equations,  existence  and  uniqueness  of  solutions  also  depends  on  the 
index  x  of  the  operator  (as  defined  in  [25]).  For  the  integral  operators  under  consideration, 
X  €  {  —  1,0, 1}.  When  x  =  the  operator  is  a  Quasi-Fredholm  operator,  and  the  Fredholm 
Alternative  is  applicable.  When  x  =  —  1,  and  the  adjoint  of  the  homogeneous  equation 
has  no  nontrivial  solutions,  then  a  solution  a  G  X^[a,c]  exist,  but  is  unique  only  up  to  a 
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constant.  When  x  =  Ij  then  the  adjoint  homogeneous  equation  has  at  least  one  nontrivial 
solution,  and  a  solution  <t  G  Z*[a,c]  exists  only  if  the  free  term  /  is  orthogonal  to  the 
solutions  of  the  adjoint  homogeneous  equation. 

When  the  operator  P  is  a  weakly  singular  first  kind  integral  operator,  then  there  may 
not  be  a  solution  c  G  L^[a,c],  since  for  this  class  of  integral  equation  P  is  a  compact 
operator.  However,  if  such  a  solution  exists,  and  the  operator  is  positive  definite,  then  the 
solution  is  also  unique  (see  [26]).  □ 

We  will  require  the  four  operators 

Paa  :  L^[a,b]-~*  L^[a,b]  , 

Pab  :  L^[b,c]-.L^[a,b], 

Pba  :  L^[a,b]^L^[b,c], 

Pbb  :  L^[b,c]-^  L^[b,c], 

defined  by 


•6 


Paa{<^)(x)  =  A -<7(1)+  /  K{x,t)  ■  cr(t)dt, 

Ja 

(3.39) 

Pab(o)(x)  =  1  K(x,t)-<T(t)dt, 

Jb 

(3.40) 

PBAi<^)(x)  =  /  K(x,t)-(T(t)dt, 

Ja 

(3.41) 

Pbb{<^){x)  =  A  •  o-(i)  +  J  K{x,t)  ■  a{t)dt. 

(3.42) 

The  integral  equation  (3.37)  can  therefore  be  rewritten  in  the  form 

Paa{<^\a)  +  P AB{o\b)  =  f\Ai 

(3.43) 

Pba{(J\a)  +  Pbb{(T\b)  =  f\B- 

(3.44) 

We  denote  subintervals  [ci,C2]  and  [di,d2]  of  [a,c]  by  C  and  D,  respectively,  and  assume 
that  C  and  D  are  disjoint.  We  define  Pcd  ■  L^{D)  — ►  L^{C)  via  the  formula 

PcDi<^){x)=  I  k{x,t)  ■  <7{t)dt,  (3.45) 

J  dl 

so  that  Pcd  is  the  operator  P  restricted  to  the  domain  C  x  D. 

Theorem  3.1  determines  the  circumstances  under  which  Pcd  may  be  represented  via 
Chebyshev  approximation.  Proofs  may  be  found  in  [21]  and  [5],  with  the  most  general  proof 
located  in  [10]. 

Theorem  3.1  Suppose  [cj,  C2],  [di, ^2]  ore  disjoint  subintervals  of[a,c],  and  let  c;  =  C2  - 
Ci.di  =  d2  —  dl,  and  r  =  min{|ci  —  d2|,|c2  —  dij}-  Suppose  further  that  Pcd  given  by 
(3.45),  with  a  kernel  of  the  form  (3.35),  and  the  interpolation  operators  are  given 
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by  (3.8),  (3.9),  respectively.  Suppose  finally  that  Chebyshev  nodes  tif,,tijy  are  defined  by 
(3.18).  If  mi, m2, given  by  Lemmas  3.3-3. 5,  then 

1.  If  Cl  <  r  (the  interval  [ci,C2]  is  well-separated  from  [di,d-^),  then 

\\PcD  -  rl’  ■  m,(PcD)ll  =  O  .  (3.46) 

2.  If  di  <  T  ( the  interval  [dj ,  ^2]  Is  well-separated  from  [ci ,  C2\),  then 

\\PCD  -  m2iPcD)  •  ^"’11  =  o  .  (3.47) 


3.  If  Cl  <  r  and  also  di  <  r  (each  of  the  intervals  [ci,C2]  and  [di,d2]  is  well-separated  from 
the  other),  then 

\\PCD  -  4’  ■  m^iPcD)  •  4>'^\\  =  o  .  (3.48) 

Let  6  =  (a  +  c)/2  denote  the  midpoint  of  the  interval  [a,c].  We  denote  the  subintervals 
[a,  6]  and  [6,c]  of  [a,  c]  by  A  and  B,  respectively. 

Remark  3.4  The  purpose  of  this  section  is  to  develop  elficient  representations  for  the 
discretized  operators  Pab,Pba-  Since  A  and  B  are  not  well-separated  intervals,  we  cannot 
apply  Theorem  3.1  to  the  operators  themselves.  However,  there  are  subintervals  of  A  and 
B  which  are  well- separated.  We  now  proceed  to  decompose  PaBjPba  into  operators  acting 
on  smaller  subintervals,  so  that  Theorem  3.1  applies  to  these  restricted  operators.  □ 

Suppose  that  we  are  given  positive  integers  and  that  r^,rg  are  defined  by  the 

equations 

rA 

tb 

If  g^i  >  1,  then  we  denote  the  subintervals  [a, 6—  r^j  and  [6  —  r^,6]  of  i4  by  Ai  and  Aq, 
respectively,  and  similarly  if  gg  >  1  we  denote  the  subintervals  [b,b+  and  [6-|-  tq,c]  of 
B  by  Bq  and  B\,  respectively.  If  g^t  >  2,  then  we  denote  g^i  -  1  subintervals  A2,.  ■ .,  Aq^  of 
>Ii,  via  the  formulae 

A.  =  [b-rA-2'-\b-ryi-2'-^}  (i  =  2, . . . ,  g^  -  1),  (3.51) 

Aq^  =  [u,6-r^-2’s-2]. 

Similarly,  if  gg  >  2,  then  we  denote  the  ge  —  1  subintervals  B2, .  ■  • ,  Bqg  of  B\  using  the 
formulae 


jb-a) 

jc-b) 

27B-1 


(3.49) 

(3.50) 


B,  =  [6+re•2'-^h-h^B•2'-']  (i  =  2. . .  .,gB  -  1). 

Bqg  =  [6+rB-2’«-2,c]. 


(3.52) 
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When  =  1,  we  define  both  Ai  and  A2  to  be  the  interval  [6,6];  when  qs  =  1,  we  define 
Bi  and  B2  to  be  the  interval  [6,6].  When  qA  =  2,  we  define  A2  to  be  equivalent  to  Ai,  and 
similarly  we  define  B2  to  be  equivalent  to  Bi  when  qs  =  2.  Thus,  for  all  values  of  qAiQB^ 

<IA 

(3-53) 

1=2 

Bx  =  (3-54) 

t=2 

Given  operators  Pab,Pba  defined  by  (3.39),  (3.40),  respectively,  we  may  express  Pab 
and  Pba  via  the  sums 

<IA  QB 

^.4b(<7|b)  =  PAoBo{<^\Bo)+^BA,Boi‘^\Bo)  +  Yl^^oBji‘^\Bj)  (3.55) 

i=2  j=2 

1A  9B 
i=2  j=2 

9a  QB 

PBAi(r\A)  =  PBoAQi(r\Ao)  +  Yl^B,Aoi(^\Ao)  +  Yl^BoA,i(^\Aj)  (3-56) 

t-2  3=2 

9b  9a 
i=2  j=2 

The  following  theorem  is  the  principal  analytical  tool  of  this  chapter.  It  permits  factor¬ 
izations  of  Pab  and  Pba  which  approximate  these  operators  to  high  accuracy.  Corollary 
3.1  then  proves  that  for  the  discretizations  of  interest,  few  coefficients  are  required  in  the 
discretization  of  these  factorizations. 


Theorem  3.2  Suppose  ^  :  L^{Ao)  x  i)  _>  L^{A)  is  given  by  the  expression 


I  ^A. 


t/’u  = 


Sa 


V^A,A-r 


Pao  / 


(3.57) 


where  Iaq  denotes  the  identity  operator  for  the  interval  Aq,  and  suppose  further  that  ^’2^  : 
L^[Bq)  X  RP  _  L^(B)  is  given  by  the  expression 
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Then 


WPAB-rPi^-MAB-tl’Ij  = 


(3.59) 


\\Pba 

-  t4’2B  •  Mba 

*U  =  o(^ 

O' 

(3.60) 

where  Mab 

:  L^{Bq)  X  RP-iiB-i)  ^  L^{Ao) 

X  Rp  («^~^)  is  given  by  the  expression 

/  rni{PA,^Bo) 

m^iPA^^B^) 

msiPA^^B,) 

"'3(Pa,^B,3 

)  >1 

mi(PA,^_,Bo) 

p^3(Pa,^_,B2 

)  rn2iPA,^_^B3) 

•  •  •  Tn2{PA^^_^B 

Mab  = 

• 

1 

’  ■ .  i 

y 

mi{PA2Bo) 

mz{PA-,B^) 

mziPA^Bz) 

•••  m3{PA2B„g) 

V  ^Aq  Bo 

^2(PaoBj) 

m2iPAoB3) 

■■■  m2{PAoB„g) 

j 

(3.61) 

and  where  Mba  ■  L‘^{Ao)  x  Rp4?.4  1)  L^(Bo)  x  Rp(9B“1)  is 

given  by  the  expression 

<  m2{PBQA„^) 

n^2(PBoA,^_i) 

m2{PBoA7) 

PBoAo  \ 

”^3(X'b2A,^) 

"»3(^’b2A,^_i) 

•••  msiPB^Ai) 

miiPsiAo) 

Mba  = 

mz{PB,A,^.l) 

mz^PB^Ai) 

mi{PB3Ac) 

y 

(3.62) 

\  Tn^{PB,^^A^^) 

m3(PB,BA,^_t 

)  rnziPB^^At 

)  rniiPB^^A^) 

and  where  mi  :  L^  x  L^  x  L^,  m2  :  L'  x  L^  L^  x  R^,  m3  :  x  RP  x  Rp  are 

defined  by  (3.24),  (3-27),  (3.30),  respectively. 


Proof.  For  each  t  =  2, . .  .,qA,  Ai  is  well-separated  from  Bq.  Therefore,  by  Theorem  3.1  we 
may  represent  operators  Pa.Bo  PboA.  to  0(1 accuracy  using  the  approximations 

^A,Bo  »  V’A,  •mi(PA.Bo). 

PboA.  «  m2iPBoA,)-i’A,^  (3-64) 

for  each  i  =  2, ...,94.  In  addition,  each  subinterval  R,  is  well-separated  from  Aq,  for 
t  =  2, . .  .,qB-  Therefore,  by  Theorem  3.1  we  may  represent  operators  PaqB,  and  Pb.Aq  to 
0(l/p*'~*)  accuracy  using  the  approximations 

PaoB.  «  »^2(PaoB.)- V’B,,  (3-65) 

Pb,Ao  ~  i’B,  ■  miiPB^Ao),  (3.66) 

for  each  i  =  2,...,qB-  For  each  i  =  2, ...,9,4  and  for  each  j  =  2,..., 9b,  each  of  the 
subintervals  Ai  and  Bj  are  well- separated  from  the  other.  Therefore,  by  Theorem  3.1  we 
may  represent  operators  Pa,Bj  and  PbjA,  toO(l/p*“')  accuracy  using  the  approximations 

Pa,B,  ~  i’A.  ■  m3{PA,Bj)  ■  i’Bj-: 

Pb,A.  ~  i'B,  •  m3{PBjA,)  ■ 


(3.67) 

(3.68) 
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for  each  i  =  2,...,qA  and  for  each  j  -  2, ...,9b.  The  approximations  (3.59),  (3.60)  are 
now  obtained  from  the  approximations  (3.63)-(3.68)  and  from  equations  (3.55),  (3.56).  □ 
The  particular  values  of  the  integers  qA  and  qg  used  in  this  section  are  discretization 
dependent.  In  this  chapter,  we  wiU  be  using  the  following  three  discretizations: 

Discretization  Cl:  Both  A  and  B  contain  n  points,  n  =  p-k,  where  A:  is  a  positive  integer, 
/lo  and  Bq  contain  p  Chebyshev  roots  respectively.  The  remaining  n  —  p  points 

for  A  and  B  are  placed  in  A2  and  B2  as  follows: 

1.  p  Chebyshev  nodes  are  contained  in  the  subinterval  [(a  +  6  —  r)/2,b  —  r]  of  A2  and  p 
nodes  are  contained  in  the  subinterval  [6  +  r,  (26  +  r)/2]. 

2.  For  z  =  2, . . . ,  A:  —  2,  both  the  subinterval  [(a  +  6  —  r)/2‘,  (a  +  6  -  r)/2*”^]  of  A2  and  the 
subinterval  [(26+  r)/2‘~^,(26+  r)/2‘]  of  B2  contain  p  Chebyshev  roots. 

3.  Finally,  the  subinterval  [a,  (a  +  6  —  r)/2*'”^]  of  A2  and  the  subinterval  [(26  +  r)/2*“^,  c] 
each  contain  p  Chebyshev  roots. 

Thi  s,  for  Discretization  Cl,  qA  =  qe  — 

Discretization  C2:  The  structure  of  /Iq  and  A2  are  the  same  as  in  Discretization  Cl. 
However,  there  are  only  p  points  in  the  interval  B,  located  at  the  p  Chebyshev  roots 
For  this  discretization,  9^  =  2,  and  qs  =  1. 

Discretization  D:  Both  A  and  B  contain  n  points,  where  n  =  p  •  2*,  for  some  integer 

k.  The  number  of  points  in  each  of  the  subintervals  and  B,  is  proportional  to  the  size 
of  the  subintervai.  Aq  and  Bq  contain  the  p  Chebyshev  roots  and  respectively, 
where  the  Chebyshev  roots  are  defined  by  (3.18).  A2  and  B2  each  contain  the  p  Chebyshev 
roots  ,  tig^,  respectively.  For  the  subintervals  A3  and  B3,  we  have  p  Chebyshev  roots  in 
each  half  of  the  subinterval,  for  a  total  of  2p  points  in  the  subinterval.  For  the  remaining 
subintervals  A,  and  Bi,  for  i  =  3,..., 9,  we  divide  each  subinterval  into  2''"^  pieces,  and 
place  the  p  Chebyshev  roots  in  each  piece  (so  that  the  total  number  of  points  is  p  •  2’~^). 
Thus,  for  Discretization  C,  94  =  9b  =  A:  +  1. 

Corollary  (3.1)  directly  follows  from  Theorem  3.2  and  from  Discretizations  Cl,  C2,  D. 

Corollary  3.1  Suppose  Mab  ond  Mba  given  by  (3.61)  and  (3.62).  Then, 

l.  When  Mab  Mba  discretized  via  Discretization  Cl,  then  the  corresponding 
discretized  operators  MaBi  Mba  require  4  •  p^  coefficients  each. 

2.  When  Mab  Mba  ure  discretized  via  Discretization  C2,  then  the  corresponding 
discretized  operators  Mab-  ^^BA  require  2  ■  p^  coefficients  each. 

3.  The  discretization  of  both  Mab  ond  Mba  cia  Discretization  D  requires  p^  +  p^  •  A:  +  p^  • 
A  +  p^  •  k^  coefficients.  Equivalently,  each  operator  Mab,  Mba  requires  p^  ■  (log2(n/p)  +1)^ 
coefficients  for  this  discretization. 


Remark  3.5  The  analytical  apparatus  of  this  chapter  strongly  resembles  the  analytical 
apparatus  of  Chapter  2.  This  is  largely  because  for  both  chapters  the  observations  of 
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Remark  2.5  apply —  that  just  as  P  was  decomposed  into  operators  PaAiPbb  low-rank 
operators  Pab^Pba,  so  we  can  decompose  each  of  Paa  and  Pbb  into  fonr  operators,  two  of 
which  are  of  low  numerical  rank.  If  this  process  is  continued,  we  obtain  the  structure  shown 
in  Figure  2.1.  However,  while  the  integral  operators  of  Chapter  2  were  of  low  analytic  rank, 
the  operators  of  the  present  chapter  are  of  low  rank  only  when  discretized.  □ 

3.4  Quadrature  Formulae  for  Singular  Integral  Operators 

In  this  section,  we  develop  quadrature  formulae  to  be  used  in  the  algorithms  of  this  chap¬ 
ter.  The  results  of  this  section  can  be  used  to  construct  px  p  discretizations  of  singular  or 
weakly  singular  integral  operators.  These  discretizations  are  well- conditioned  (the  under¬ 
lying  quadrature  weights  are  nearly  all  positive),  and  exhibit  pth  order  convergence  to  the 
analytic  integral  operator. 

3.4.1  Product  Integration  for  Singular  Functions 

We  first  consider  the  numerical  evaluation  of  integral  equations  of  the  form  (3.1),  with  the 
kernel  k  given  by  (3.35),  and  the  singularity  s  given  by  either  (3.2)  or  (3.3). 

Remark  3.6  The  method  of  this  section  is  not  required  for  the  Cauchy  singularity  (3.4),  as 
high-order  quadrature  formulas  for  this  singularity  already  exist  (see,  for  example,  [8], [9]). 
□ 


We  denote  an  n-point  discretization  of  the  interval  [—1,1]  by  [xj,  X2, .  ■ . ,  x„].  Given  this 
discretization  and  quadrature  weights  Wij  (1  <  i,jf  <  n),  the  Nystrom  algorithm  replaces 
the  integral  equation  (3.1)  with  the  linear  system 

n 

A  ■  <7,  -f  ^  w,j  •  {ki{xi,Xj)s{xj  -  X.)  A  k2{xi,Xj))  ■  =  f{xi)  (1  <  t  <  n).  (3.69) 

j=i 

The  solutions  <ti,  . . . ,  cr„  are  then  viewed  as  approximations  to  the  solution  a  of  (3.1)  at  the 
nodes  xi,...,x„.  As  is  well-known,  if  (3.1)  has  a  unique  solution,  then  for  a  wide  class  of 
quadrature  formulae  and  sufficiently  large  n,  the  system  (3.69)  has  a  unique  solution,  and 
the  solution  of  (3.1)  converges  to  the  analytic  solution  at  the  same  rate  as  the  underlying 
quadrature  formula. 

Unfortunately,  standard  quadrature  formulae  (for  example,  end-point  corrected  trape¬ 
zoidal  quadrature  formulae,  or  Chebyshev  quadrature)  exhibit  extremely  poor  convergence 
for  equations  of  the  form  (3.1),  due  to  the  singularity  s.  Appropriate  quadrature  formulae 
are  usually  obtained  by  some  form  of  product  integration  (see,  for  example,  [17]).  Such 
formulae  were  obtained  in  [27],  [3]  for  end-point  corrected  trapezoidal  quadrature  rules;  we 
now  use  the  method  of  these  papers  to  construct  weights  for  Chebyshev  quadrature. 

We  assume  that  (3.1)  has  been  discretized  at  the  2p  zeroes  [xi,...,X2p]  of  the  2pth 
Chebyshev  polynomial  T^p.  Using  the  notation  of  [32],  we  define  generalized  moments  of 


56 


CHAPTER  3.  ONE-DIMENSIONAL  INTEGRAL  EQUATIONS 


the  moment  function  s,  with  s  singular  at  Xi,  by  the  expression 


=  J  ^  7j-i(0  ■  •*{<  -  **)  dt  (1  <  j  <  p). 


(3.70) 


For  each  i,  we  now  consider  the  following  system  of  linear  algebraic  equations  with  respect 
to  the  unknowns  Wij  {1  <  i  <  n): 


for  j  =  1 , . . . ,  p,  and 


2P  fl 

V  Tj-i(xk)  ■Wik=  /  Ij-iCO  dt, 
•'-1 


(3.71) 


^  ^  Tj~p—\(^X k^  ■  5(3?^  Xi')  ‘  Wik  —  TTlj—p(^Xi)y 


(3.72) 


for  j  =  p  +  1 , . . . ,  2p.  Obtaining  the  quadrature  weights  Wij  therefore  requires  the  solution 
of  a  total  of  2p  linear  systems. 

Remark  3.7  While  the  2p^  moments  7nj(a:,)  can  be  determined  analytically,  this  is  quite 
burdensome  for  large  p.  Instead,  we  evaluate  the  moments  numerically,  using  adaptive 
gaussian  quadrature  (see,  for  example,  [30]).  The  expression  (3.70)  can  be  evaluated  to 
high  accuracy,  provided  the  singularity  is  given  by  (3.2),  or  is  given  by  (3.3)  with  a  positive 
exponent. 

If  the  singularity  is  given  by  (3.3)  with  negative  exponent,  then  adaptive  gaussian 
quadrature  will  yield  high  accuracy  only  if  the  location  of  the  singularity  is  zero  (for  other 
locations  of  the  singularity,  the  subinterval  divisions  can  not  be  made  sufficiently  fine  to 
yield  high  accuracy).  However,  the  following  expression  is  equivalent  to  (3.70),  and  can  be 
evaluated  accurately  using  adaptive  Gaussian  quadrature: 


ij{xi)=f  Tj-i{xi-t)-s(t)dt  (l<j<p). 
Jx,-1 


(3.73) 


The  existence  of  unique  solutions  to  the  linear  systems  and  the  convergence  rate  of  the 
quadrature  formulae  are  given  by  the  following  two  lemmas,  which  were  presented  in  [3]  (in 
a  slightly  different  form). 

Lemma  3.6  For  each  Xi  {I  <  i  <  2p),  the  linear  system  (3.71)-(3.72)  has  a  unique  solu¬ 
tion. 

Lemma  3.7  For  each  x,,  the  convergence  rate  of  the  quadrature  weights  Wij  {I  <  j  <  2p) 
is  at  least  p. 

Unfortunately,  as  the  order  of  the  quadrature  formulae  p  increases,  the  linear  systems 
given  by  (3.71)-(3.72)  become  increasingly  ill-conditioned.  As  a  result,  the  weights  Wij 
become  large,  to  the  point  that  substantial  accuracy  is  lost  when  p  >  8. 
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3.4.2  High-Order  Quadrature  Formulae  with  Well-Conditioned  Weights 

The  algorithms  of  [27],  [3]  have  favorable  analytic  properties,  but  produce  unacceptably 
large  quadrature  weights  for  an  order  of  convergence  p  >  8.  In  contrast,  the  algorithm  we 
now  present  is  a  series  of  linear  algebraic  fixes  to  [27],  [3].  We  have  no  proofs  concerning  the 
performance  of  the  algorithm  in  this  section,  but  in  practice  it  has  been  used  to  generate 
extremely  high-order  methods  with  quadrature  weights  that  are  nearly  aU  positive. 

We  rewrite  the  system  (3.71)-(3.72)  in  the  form 

Miwi  =  /i,  (3.74) 

where  Mi  €  is  the  system  of  Chebyshev  polynomials  and  Chebyshev  polynomials 

times  a  singular  function,  wi  €  is  the  vector  of  quadrature  weights,  and  /  6  R^**  is 
the  vector  of  analytic  integrals.  Assume  that  to  some  precision  c,  the  rank  of  the  matrix 
Ml  :  f  (the  matrix  Mi  with  an  additional  column  /)  is  n.  Then,  an  orthonormal  matrix 
Qi  €  R2p><2p  jjg  constructed  which  maps  the  left  nuUspace  to  zero  (to  precision  c). 
More  precisely,  we  obtain  an  orthonormal  matrix  Qi  such  that  rows  n  -f-  1, . .  .,2p  of  the 
matrix  Q  •  Mi  are  zero,  to  precision  c.  Let  Q2  €  R”^^p  denote  the  matrix  with  its  n  rows 
equivalent  to  the  first  n  rows  of  Qi.  Then,  to  precision  c,  the  linear  system 

Q2M1W1  =  Q2/1  (3.75) 

is  equivalent  to  the  linear  system  (3.74).  For  convenience,  we  rewrite  the  system  (3.75)  as 
the  system 

M2W1  =  f2,  (3.76) 

with  M2  €  R'^’^^P  defined  by 

M2  =  <52- Ml,  (3.77) 

and  /2  €  R"  defined  by 

f2  =  Q2-fi.  (3.78) 

By  Lemma  3.6,  there  exist  solutions  to  (3.76).  We  obtain  Wi  by  solving  the  following 
least  squares  problem  with  linear  constraints  (see,  for  example,  [15]): 

Minimize  |itni||2  subject  to  M2'Wi  =  /2.  (3.79) 

While  solving  the  system  (3.79)  does  provide  solutions  for  p  >  8,  the  weights  wi  obtained 
are  stiU  quite  large.  However,  if  we  construct  weights  for  points  located  at  the  zeroes  of 
a  Chebyshev  polynomial  n  >  2p,  these  weights  can  be  nearly  all  positive.  The  original 
equations  (3.71)-(3.72)  thus  become,  for  each  xi: 

r,_i(x*)  ■  w,,  =  /'  Tj.i{t)dt,  (3.80) 

k=i 

for  j  =  1. ....  p,  and 

n 

Y^Tj-p-i{xk)  ■  s{xk  -  X,)  ■  Wik  =  mj_p(x,), 
k-i 


(3.81) 
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foT  j  =  p  4-  1, . . 2p.  Our  experiments  indicate  that  when  n  >  3p,  the  weights  are  either 
all  positive  or  nearly  all  positive.  Unfortunately,  we  have  no  explanation  for  this  behavior. 

Remark  3.8  In  [3],  Alpert  interpolates  a  function  discretized  at  equispaced  nodes  to  a 
Chebyshev  discretization,  in  order  to  delay  the  growth  in  weights.  A  similar  interpola¬ 
tion  strategy  can  be  used  here  to  generate  a.  p  x  p  integral  operator  with  pth  order  con¬ 
vergence,  even  when  the  number  of  quadrature  weights  required  for  a  particular  row  is 
large  (we  use  3p  weights  per  row).  Assume  an  operator  is  to  be  applied  to  a  function 
/  =  [/(^i[a,ci’/(^a,ci)»-  ^^^h  givcu  by  (3.18),  and  assume  further¬ 

more  that  /  6  C^[a,c].  Denote  by  M  the  p  x  3p  matrix  containing  the  discretized  kernel 
function  cl?  J/j)  ®3.ch  ^j,  the  points  y\,y2,  •  •  •,  l/3p  are  chosen  in  some  convenient 

fashion).  Let  C  :  R'’  — *  Rp  denote  the  discrete  Chebyshev  transform  operator  of  dimension 
p.  Then,  M  •  f  will  yield  the  p  Chebyshev  coefficients  of  the  function  /.  These  coefficients 
can  then  be  used  to  yield  the  interpolated  values  of  /  at  the  points  yi,  j/27  •  •  •  >  ysp- 

By  itself,  this  interpolation  procedure  is  of  little  interest.  However,  when  combined  with 
the  algorithms  of  this  chapter,  the  order  of  convergence  of  the  discretized  integral  operator 
is  improved  from  p/3  to  p.  □ 

3.5  The  Analytical  Apparatus  for  Singular  Solutions 

In  the  following  two  sections,  we  use  the  results  of  Section  3.3  to  produce  the  apparatus  for 
the  rapid  solution  of  the  integral  equation  (3.1).  By  Corollary  3.1,  the  exact  representation 
of  the  operators  Pab  and  Pba  defined  by  (3.40),  (3.41)  is  discretization  dependent.  In  this 
section,  we  present  results  for  Discretizations  Cl  and  C2 —  these  discretizations  are  used 
when  there  is  an  end-point  singularity  in  the  solution  of  (3.1).  In  Section  (3.6),  we  present 
results  for  Discretization  D;  this  discretization  is  used  when  the  solution  to  (3.1)  is  smooth. 

The  fast  algorithms  of  this  chapter  are  based  on  the  fund2unental  observation  that  the 
solution  to  the  integral  equation  (3.1)  on  the  entire  domain  [o,c]  can  easily  be  constructed 
from  the  solution  of  the  two  independent  integral  equations  (3.43)-(3.44).  one  defined  on 
[a,  6]  and  one  on  [6,  c].  This  leads  naturally  to  a  recursive  algorithm,  in  which  independent 
solutions  on  a  large  number  of  subintervals  are  successively  merged  until  the  full  solution 
is  obtained.  A  precise  formulation  of  the  construction  and  the  resulting  numerical  scheme 
will  require  some  notation. 

3.5.1  Notation 

We  let  b  =  (a  +  c)/2  denote  the  midpoint  of  [a,c],  and  denote  the  subintervals  [a,  6]  and 
[5.  c]  by  A  and  B,  respectively.  In  addition,  we  require  the  operators 

P  :  I^[a,c]  —  T^[a,c], 

Paa  ■  L^[a,b]-^  L'^[a,b], 

Pab  ■  L^[b,c]-^  L'^[a,b], 

Pba  :  L^[a,b]^  L^{b,c], 


3.5.  THE  ANALYTICAL  APPARATUS  FOR  SINGULAR  SOLUTIONS  59 

Pbb  :  L^[b,c]^L^[b,c], 

defined  by  (3.33),  (3.39),  (3.40),  (3.41),  (3.42),  respectively.  Due  to  Theorem  (3.2),  if  we 
are  using  Discretization  Cl  then  we  may  approximate  Pab  and  Pba  via  the  expressions 


Pab  «  •  M\b  ■  V’ls, 

Pba  «  Vib  ■  ^BA  ■ 


(3.82) 

(3.83) 


:  W’ X  L\Ao) L\A), 

X  L\Bo)  ^  L\B), 

the  interpolation  operators  of  the  form  (3.57),  given  by  the  expressions 


II 

^  i’A2 

Iao 

i’lB  = 

/  V’B2 

IBo 

(3.84) 

(3.85) 


given  via  the  formulae 


M\b  :  HTxL\Bo)^BTxL'^{Aol 

M^ba  ■  X  L\Ao) R”  X  L\Bo), 


^BA  -  ^ 


miiPA^Bo) 

^siPAiBi) 

(3.86) 

PaqBq 

m2{PAoB3)  )  ’ 

‘m2{PBoA2) 

PBoAo  J 

(3.87) 

m3{PB2A2) 

mi{PB2Ao)  )  ' 

Similarly,  if  we  are  using  Discretization  C2  then  using  Theorem  (3.2)  we  may  approximate 
Pab  and  Pba  via  the  expressions 


V'u  -^AB^ 


(3.88) 

(3.89) 


MIb  ■  LHB)^R^xL\Ao), 
MIa  ■  R’’ X  L\Ao)  ^  L^iB), 
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given  via  the  formulae 


mIb 

_  /  rni{PA^B)  \ 

Pa,Bo  )  ’ 

(3.90) 

MIa 

=  (  Tn2{PBA2)  PbAo  )  • 

(3.91) 

We  define  the  operator 

Q  :  L\a,c]  x  Rp  L\a,c\  x  RP 

by  the  expression 

m/' 

Qix){x)  =  xix)+  k{x,t)-x{t)dt. 

J  a 

(3.92) 

W'e  additionally  require  the  four  operators 

Qaa  '■ 

l2[a,6]xRP-i2|^,6]xRP  , 

Qab  '■ 

L\b,c\  X  RP  ^  L\a,b]  X  RP  , 

Qba  ■ 

L^[a,b]  X  RP  ^  T^[6,c]  x  RP  , 

Qbb  '■ 

I2[6,c]xRP->T2[6,c]xRP  , 

defined  by  the  formulae 

<3.4/t(x)(a; 

)  =  X(2^)+  /  k{x,t)-x{t)dt, 

J  a 

(3.93) 

)  =  k{x,t)-x{t)dt, 

(3.94) 

<5ba(x)(3; 

)  =  /  k(x,i)-xit)dt, 

Ja 

(3.95) 

Qbb{x){^ 

)  =  X(3-)  +  ^  k{x,i)  ■  x{t)dt. 

(3.96) 

Given  a  function  /  €  L^[a,c],  we  will  follow  the  convention  of  denoting  its  restriction  to 
.4  and  B  by  /j^  and  /|g,  respectively.  Similarly,  given  a  function  £  L^[a,c]  x  R’’,  we  will 
denote  its  restriction  to  A  and  B  by  4'\a  V’|b,  respectively. 

Given  an  interval  [61,62]  C  R  and  operators  ;  R"  — ►  L^[bi,b2],  X  •  R"*  •^>^[61, 62], 

and  let  t/’,(x)  denote  the  ith  component  of  and  let  denote  the  yth  component 

of  Then  we  define  the  inner  product  a  €  LCR”*^")  by  the  expression 

(o),_,  =  ^  'tl’dt)  ■  X]{t)  dt,  (3.97) 

where  (a),^  denotes  the  entry  in  the  ith  row  and  jth  column  of  a.  Similarly,  we  define  the 
transpose  of  c  by  the  formula 

uAt)Xjit)dt>  (3.98) 

Jb] 
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Given  interpolation  operators 


i’la.c] 

W>-^L^[a,c], 

0/i 

RP  L^{A), 

0B 

RP  X^(R), 

defined  by  (3.14),  we  will  in  this  section  refer  to  them  by  •03,  03^,  03^,  respectively. 

Given  Discretizations  Cl  and  C2,  we  require  the  interpolation  operator 

0i;R'’L"(B)->L''[a,c], 


given  by  the  formula 


We  in  addition  require  the  zero  operators 

0^  :  L^[a,c]-^L\A), 
Ob  :  RP-*i2(B). 


(3.99) 


(3.100; 


The  operators  0i  and  03  are  related  to  03^, 03^  via  the  expressions 

01m  ^  0^  ’  (3.101) 

0i|B  =  Ob  X  /b  ,  (3.102) 

03|^  =  03m  •  Cd  ,  (3.103) 

V’3,b  =  03b  •  Cu  ,  (3.104) 

where  Co,  Cu  •  R'’  — *  R^  denote  the  shifting  matrices  defined  by  (3.21),  (3.22),  respectively. 
Assuming  that  the  operators  Paai  Pbb  are  nonsingular,  we  define  the  functions 

rjA  ■■  A  ^  R, 

77b  R  — ‘  R, 

via  the  formulae 

riA  =  Pa\U\a).  (3.105) 

r}B  =  /’bM/ib).  (3.106) 

Similarly,  assuming  that  the  operators  P,  Paa-  Pbb  are  nonsingular,  we  then  define  the 
operators 

XI  :  X  L^{B) L'^[a,c]  , 

X3  ;  R’’  — ►  L^[a,  c]  , 

:  RP  X  Z‘(.4o)  -  Z"(A)  , 
d>3,  :  RP  ^  L^i.A)  , 

d3«  :  K^-L^B), 


62 


CHAPTER  3.  ONE-DIMENSIONAL  INTEGRAL  EQUATIONS 


via  the  formulae 


xi  =  P  ‘(V’l) , 

(3.107) 

X3  =  P-\^3), 

(3.108) 

=  PaaI^Ia)  . 

(3.109) 

^A  =  Paa{‘^^a)  ’ 

(3.110) 

—  P Bfi(^3s)  • 

(3.111) 

Finally,  given  the  transpose  defined  by  (3.98),  we  define  operators 


by  the  formulae 


and  the  functions 


by  the  expressions 


On 

R'’ 

X  L'^{Aq)  — »•  RP  X 

L^Ao) 

<^13 

RP 

^  RP  X  L'^{A(,)  , 

<^31 

RP 

X  L'^iAo)  -  RP  , 

RP 

-RP  , 

"11 

RP 

X 

t 

M 

X 

L^B)  , 

013 

RP 

-  RP  X  L‘^{B)  , 

«31 

RP 

X  X2(B)  _  RP  , 

O33 

RP 

-R" , 

Oil  - 

’  «13  = 

■^A  > 

031  = 

4>Ia  »  =  V’3^ 

•  4>33  . 

an 

- 

•  Xi  ^  ai3  =  i’f  • 

X3  , 

a3i 

= 

•  Xl  ,  033  =  ^3  ■ 

X3  , 

St 

€  W’xL\Ao) 

6  RP, 

€  RP  X  L\B)  , 

h 

€  RP  , 

6-^  =  ■  Va  .  ^3  =  ^3^  •  VA  , 


(3.112) 


(3.113) 


^1  =  pf  ■  , 

where  a  is  the  solution  to  equation  (3.37). 


S3  =  ibj  ■  a  , 


(3.114) 
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3.5.2  Recursive  solution  of  the  integral  equation 

We  now  consider  the  original  integral  equation  (3.37) 

P(j  =  f. 

The  main  results  of  this  section  are  the  following  two  lemmas,  which  construct  the 
solutions  a  of  equation  (3.37)  from  the  solutions  t)a,t)b  of  equations  (3.105)  and  (3.106). 
Since  the  proofs  are  quite  similar,  we  only  present  the  proof  for  Lemma  3.8. 

Lemma  3.8  If  Discretization  Cl  is  used,  so  that  operators  PaBiPba  cun  be  approximated 
via  (3.88)-(3.89),  and  if,  in  the  notation  of  Section  3.5.1,  all  six  operators  P,  Paa,  Pbb, 
Qt  Qaa>  Qbb  non-singular,  and  if  furthermore  the  operator 

Ai  ;  RP  X  L\Ao)  Rp  x  L\Ao) 

given  by 

^1  =  ^RPxL^(Ao)  ~  ■  ^AB  •  “fl  ■  ^BA‘> 

is  also  nonsingular,  then  to  accuracy  0(l/p^~^), 

^{A  =  VaF  <P1a  ■  ■  (afi  •  -^i^- 

=  VB  A  4>Ib  ■  Mb  A' ' '^11' ^AB  ' 

{bf -af,.  Ml  6^). 

Proof  Using  definitions  (3.39)-(3.42),  the  integral  equation 

Pa  =  f 

can  be  rewritten  in  the  form  (3.43)-(3.44).  The  approximations  (3.82)  and  (3.83)  for  Pab 
and  PbAi  respectively,  can  then  be  used  to  obtain  an  explicit  solution  to  the  coupled  equa¬ 
tions  (3.43)  and  (3.44)  in  terms  of  the  functions  defined  by  (3.105),  (3.106), 

(3.109),  (3.110),  respectively.  Indeed,  applying  the  operator  P^l  to  equation  (3.43)  and  the 
operator  Pgl  to  equation  (3.44),  we  have 

a|,4  +  PJ]  -PABi^^lB)  =  PlAiflAh  (3-118) 

PEh  ■  PbaIc^ia)  +  <^\B  =  PEhiflB)-  (3-119) 

Substituting  the  approximations  (3.82)  and  (3.83)  into  (3.118)  and  (3.119)  yields  the  for¬ 
mulae 

^lA  +  PEa  -i'tA-  ^^AB  ■  ■  ^\B  = 


(3.115) 

(3.116) 

(3.117) 


(3.120) 
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Pbb  ■  •  ^BA  ■  ■  (^\A  +  <7|b  =  rjB,  (3.121) 

or 

<^\A  +  <t>iy^- =  Va,  (3.122) 

<t>iB  •  ^BA-i’u-(^\A+<^\B  =  VB,  (3.123) 

where  we  have  used  the  definitions  (3.109),  (3.110)  for  d>i^  aJid  respectively.  Now, 
multiplying  (3.123)  by  <t>\^  •  M\q  ■  and  subtracting  it  from  (3.122),  we  obtain 

ih^A)-4>lA- ^AB-'<i’lB-^B-^BA-‘>l^\J-(^\A  =  (3.124) 

VA  -  4>\a  •  ^AB  •  •  VB- 

Similarly,  multiplying  (3.122)  by  d>ig  •  Mba  ‘  «^Dd  subtracting  it  from  (3.123)  results  in 
the  equation 

{Il‘^(B)-4>\b-^ba-'^\a-<^a-^ab-'1’\b)'^\b  =  (3.125) 

t]b  -  4>Ib  ■  M^a  ■  'PIa  ■ 

Due  to  (3.112)  and  (3.180),  we  can  rewrite  these  equations  in  the  form 

ih^A)  -  4>1a  ■  ^AB-o^n  ■  M^A'i^lA))'<^\A  =  VA-  4>1a  '  ^AB  '  ^  (3.126) 

ihHB)  -4nB-MlA-  («A  •  M\b  ■  i’Js))  •  =  ^B-<hB-  MhA  ■  (3-127) 

By  application  of  Lemma  2.5,  we  obtain 

<^1-4  =  ih^^A)-^  <hA- ^AB'<^\\' ^^BA-  (3.128) 

{IrpxLHAo)  ~  "A  •  M\b  •  afi  •  •  ilA  -  4>Ia  ■  M\b  •  <5f ), 

^\B  =  {h'^(B)E(t>\B- Mba'{I'Rj>xL‘^(Ao)~^\\'^AB'^\\'^Ba)~^  ‘  (3.129) 

On  •  ^^AB  ■  ^1b)  ■  ('IB  “  •  ^BA  •  ^l*)- 

The  equations  (3.116),  (3.117)  are  now  obtained  from  equations  (3.128),  (3.129)  and  equa¬ 
tion  (3.115).  □ 


Lemma  3.9  If  Discretization  C2  is  used,  so  that  operators  PaBi  Pba  be  approximated 
via  (3.82)-(3  83),  and  if,  in  the  notation  of  Section  3.5.1,  all  six  operators  P,  Paa,  Pbb, 
Q-  Q  a.a<  Qbb  ore  non-singular,  and  if  furthermore  the  operator 

A2  :  L‘^{B)  -  L'^{B) 


given  by 

A2  =  Ib  -  Pbb  '  ^^BA  ■  ®n  ■  ^\b'  (3.130) 

is  also  nonsingular,  with  defined  by  (3.90)-(3.91),  then  to  accuracy  0(11 

<^\A  =  Va  4>Ia  '  ^^AB  '  i^2  *  ■  Pbb  '  ^^ba  '  (3.131) 

(<5i*-On  ■  MIb  ■  Vb)  -  Vb), 

(7|fl  =  2i^^iVB~PBh-MBA-bt)-  (3.132) 
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3.5.3  Further  Analytical  Results  for  Discretization  C2 

We  now  collect  a  number  of  identities  which  are  necessary  for  the  algorithm  to  be  presented 
in  Section  3.7.  First,  we  apply  Lemma  3.9  to  the  particular  cases  /  =  V’l,  /  =  V’a  to  obtain 
analytical  expressions  for  the  functions  xi  and  X3  defined  in  equations  (3.107)  and  (3.108). 
We  omit  the  proofs  here,  as  they  are  quite  similar  to  the  proofs  for  Lemmas  3.8-3. 9 

Corollary  3.2  If  Discretization  C2  is  used,  so  that  operators  PaBiPba  ^  approxi¬ 
mated  via  (3.82)-(3.83),  and  if,  in  the  notation  of  Section  3.5.1,  all  six  operators  P,  Paa, 
Pbb,  Q,  Qaa>  Qbb  are  non- singular,  then  to  accuracy  0{\lp^~^), 


Xl|>  —  '  ^AB  '  ^2^  ’  ^BB  '  ^BA' °‘13)  (3.133) 

(-<^1^  •  •  (Aj ^  •  Pgg  ■  ■  a^j  •  M\q  ■  Pgg  +  Pbb))  ’ 

=  (-A2-'-Pi^-M|^.a^3)x(A2-'.PB^),  (3.134) 

X3|^  =  ■  Cd  +  ■  NI^B  ‘  {^2^  '  ^BB  ‘  ^BA  '  (3.135) 

(^13  •  -  «n  •  ^AB  •  <hB  •  (^u)  -  4>3b  ■  Cu)  , 

X3,«  =  ^^^■{(hB-Cu-PBB-MlA-O^X3-CD).  (3.136) 


where  the  coefficients  af^  and  are  given  by  equation  (3.112). 

We  will  also  require  analytical  expressions  for  the  inner  products  ,  and  8z  defined  in 
(3.114)  in  terms  of  the  restricted  inner  products  and  defined  in  (3.114).  The  proofs 
follow  readily  from  Lemma  3.9. 

Corollary  3.3  If  Discretization  C2  is  used,  so  that  operators  PaBiPba  ^oy  be  approxi¬ 
mated  via  (3.82)-(3.83),  and  if,  in  the  notation  of  Section  3.5.1,  all  six  operators  P,  Paa, 
Pbb,  Q,  Qaa,  Qbb  ore  non-singular,  then  to  accuracy  0{\/p'^~^), 

h  =  if  •  =  (V3^  X  Ofl)^  •  <T|4  +  (O4  X /b)^  •  cr|B  (3.137) 

^3  +  ®31  ■  ^^AB  ■  (^^2  ^  ■  ^BB  '  ^BA  '  ~  “n  ■  ^AB  ' ~  ^b) 

^\B 

h  =  Clr-rl^^^-CT\B  +  Cl-8^  +  Cl-ai,-MlB-  (3-138) 

■  Psh  ■  -  On  •  Mis  ■  riB)  -  Vb)  • 

Special  cases  of  Corollary  3.2  are  obtained  when  /  =  or  /  =  1P3.  The  proofs  follow 
easily  from  the  definitions  of  \i  and  \3  in  (3.107)  and  (3.108). 
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Corollary  3.4  If  Discretization  C2  is  used,  so  that  operators  Pab^Pba  ^  approxi¬ 
mated  via  (3.82)-(3.83),  and  if,  in  the  notation  of  Section  3.5.1,  all  six  operators  P,  Paa, 
Pbb,  Q,  Qaa>  Qbb  non-singular,  then  to  accuracy  0{ljp'‘~^), 


,  T  (  Mpp  X  Npj)  \ 

on  j. 

with  .Mpp,  Npp  ;  given  by 

Mpp  =  +  031  ■  '  ^2^ '  Pbb  ’  ^^ba  '  “is  > 

Npp  =  -c^i:Mls-i^2'PEh-MlA-c^uMAB-PEB  +  PBB)- 


(3.139) 


Corollary  3.5  If  Discretization  C2  is  used,  so  that  operators  Pab^Pba  be  approxi¬ 
mated  via  (3.82)-(3.83),  and  if,  in  the  notation  of  Section  3.5.1,  all  six  operators  P,  Paa, 
Pbb-  Q-  Qaa,  Qbb  ore  non-singular,  then  to  accuracy  0(1  fp'‘~^), 

ai3  =  V-r  ■  X3  =  ^  j  ,  (3.140) 


with  .Mpp  :  R^  — ‘  RP  given  by 

Mpp  =  o^  •  Cd  +  031  •  ■  (Aj  ^  •  Pbb'  ^ba  ' 

(<>13  •  Qd  -  On  •  Mab  ■  <hs  ■  Cu)  -  4>Zb  ‘  ^v)  ■ 


Corollary  3.6  If  Discretization  C2  is  used,  so  that  operators  Pab,Pba  ^oy  be  approxi¬ 
mated  via  (3.82)-(3.83).  and  if,  in  the  notation  of  Section  3.5.1,  all  six  operators  P,  Paa, 
Pbb-  Q,  Qaa,  Qbb  ore  non-singular,  then  to  accuracy  0(1  fp^  *), 

031  =  v^-\x=  (3.141) 

■  '^'Ib  '  4-  Cq  ■  031 '  ^AB  ’  ^2^  '  Pbb  '  ^^ba  '  “fa)  ^ 

(-Cq  -  031  .  M^b  •  {‘^2^  Pbb'  ^^ba  '  On  •  ^ab  Pbb^  Pbb))  • 


Corollary  3.7  If  Discretization  C2  is  used,  so  that  operators  Pab,Pba  ^oy  be  approx¬ 
imated  via  (3.82)-(3.83),  and  if.  in  the  notation  of  Section  3.5.1,  all  six  operators  P, 
Paa-Pbb-  Q-  Qaa-  Qbb  ore  non-singular,  then  to  accuracy  0(1/ p^'~^), 

033  =  ^'3  ■  X3  =  (3.142) 

("r  ■  ^'3b  ■  X3|b  +  Cp  ■  o^  -  Cd  +  Cp  ■  031  .  •  (A2  ^  •  PgB  '  ^^ba  ' 

(o'u  ■  Cp  -  ofi  •  M\b  ■  03g  Cl  )  -  03g  •  Cv). 


Finally,  combining  Lemma  3.9  with  the  expressions  (3.133)-(3.136),  w'e  have 
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Corollary  3.8  Suppose  Discretization  C2  is  used,  so  that  operators  Pab,Pba  he  ap¬ 
proximated  via  (3.82)-(3.8S),  and  suppose  that  in  the  notation  of  Section  3.5.1,  all  six 
operators  P,  Paai  Pbb,  Q,  Qaa,  Qbb  ore  non-singular.  Suppose  further  that  the  function 
F  :  [a,c]  ^  R  is  defined  by  the  formula 

^  +  X3(a;)  •  A3  +  CT(a:),  (3.143) 

with  X2  G  and  Ai,  A3  €  R^.  Then  on  the  interval  [a, 6],  to  accuracy  0(l/p*~^), 

Fix)  =  <f>iy^ix)  ■  Pi  +  <hA^)  •  M3  +  T}Aix)  ,  (3.144) 

with  the  functions  pi  6  R^  x  L^iAo),  p^  €  R''  defined  by  the  formulae 

Pi  =  •  (Aj  *  •  Pgg  •  •  (0^*3  •  (Ai  +  C£)  •  A3) - 

0‘ll  •  ■  (PgQ  ■  A3  +  </)3g  •  Cu  ■  A3  +  T)b)  +  S^)  — 

FbB  ■  A2  -  (he  -  Cu  ■  X3-  T)b)  , 

P3  =  Ai  +  Cb  ■  A3.  (3.145) 

Similarly,  on  the  interval  [i,c],  to  accuracy  0(l/p*“^), 

Fix)  =  Xi|B(a;)  •  (  A2  )  ■'■  A3|b(^)  •  ^^3  +  (x^Bix).  (3.146) 

Proof.  Restricting  (3.143)  on  the  subintervaJs  A,B  of  [a,c],  respectively,  we  have 

F\a  =  Xi|^  *  ^  ^  +  X3|x  •  A3  +  ^  (3.147) 

F\b  =  XiiB  •  ^  j  +  X3|b  •  A3  +  (7|b  .  (3.148) 

Now,  (3.148)  is  equivalent  to  (3.146),  and  the  expression  (3.145)  results  from  combining 
(3.147)  with  (3.131),  (3.133),  (3.135),  and  comparing  the  resulting  expression  with  the 
expressions  (3.144).  □ 


Fix)  =  xi(i) 


3.6  The  Analytical  Apparatus  for  Smooth  Solutions 

In  this  section,  we  present  results  for  the  solution  of  the  integral  equation  (3.1),  when  the 
solution  is  smooth.  Like  Section  3.5,  the  techniques  of  the  present  section  are  used  to  merge 
independent  solutions  on  a  large  number  of  subintervals  until  the  full  solution  is  obtained. 
In  addition,  the  notation  and  theorems  of  the  present  section  strongly  resemble  those  of 
Section  3.5,  but  with  two  important  differences.  First,  the  underlying  discretizations  are 
different  (Discretization  D  for  the  present  section,  Discretizations  Cl  and  C2  for  Section 
3.5),  resulting  in  different  factorizations  of  Pab  and  Pba-  Second,  because  the  two  sections 
assume  different  merging  strategies,  the  relations  (3.162)-(3.165)  in  the  present  section  are 
different  from  the  relations  (3. 101)-(3. 104)  of  Section  3.5. 
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3.6.1  Notation 

We  let  6  =  (a  +  c)/2  denote  the  midpoint  of  [a,c],  and  denote  the  subintervals  [a,  6]  and 
[b,c]  by  A  and  B,  respectively.  Given  the  integer  q  associated  with  Discretization  D  (see 
Corollary  3.1),  we  define  r  =  (c  —  a)f2^,  and  define  subintervals 


Al 

=  [a,a  +  r], 

(3.149) 

Ac 

II 

1 

(3.150) 

Be 

=  [6,6+ r]. 

(3.151) 

Bf 

=  [c-r,  c]. 

(3.152) 

In  addition,  we  require  the  operators 


P 

L\a,c]  — 

L\a,c\, 

Pa  A 

L'^[a,b]-^ 

L^[a,b], 

Pab 

L\b,c]-* 

L^[a,b], 

Pba 

L\a.b]-^ 

i2[6,c]. 

Pbb 

L2[6,c]- 

I2[6,c], 

defined  by  (3.36),  (3.39),  (3.40),  (3.41),  (3.42),  respectively.  Due  to  Theorem  (3.2),  when 
we  are  using  Discretization  D  we  may  approximate  Pab  aJtd  Pba  via  the  expressions 

Pab  «  V’u  •  Mab  •  (3.153) 

Pba  «  V2b-  MBA-i’x^',  (3.154) 

with 

X  L\A,)  ^  L\A), 

X  1^(5/)  -  1^(5), 

^2^  :  L\A), 

i^2B  ■  l'(Be)  X  -  1^(5), 

interpolation  operators  given  by  (3.57)  and  (3.58),  and  with 

Mab  ■■  L'^{Bc)  x  RP<<?-i)  _  rp(?-i)  x  L'^{A^), 

Mba  ■  X  l2(/lc)  -  i''(5c)  X 

defined  by  (3.61),  (3.62),  respectively. 

We  define  the  operator 

Q  ;  L'^\a,c]  x  RP’^  -  L\a,c\  x  RP’, 

Q(\)(2')  =  \(i)  +  /  k(xj)-x{t)dt. 

J  a 


by  the  expression 


(3.155) 
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We  additionally  require  the  four  operators 


Qaa 

L^[a,b]  X  rH9-i)  _  i2^a,b]  x  Rp(’-i)  , 

Qab 

L%c]  X  RP  (9-i)  ^  ^  RP  (9-i)  ^ 

Qba 

L^[a,b]  X  Rp(9-D  _  L^[b,c]  X  RP  (9-i)  , 

Qbb 

L'^[b,c]  X  rH<j-i)  i-i[b,c]  X  RP-(«-^)  , 

defined  by  the  formulae 

r 

Qaa{x){x)  =  >:(ar)+  /  k{x,t)  ■  x{t)  dt, 

Ja 

(3.156) 

Qab{x){i)  =  ^  k{x,i)-x{t)dt. 

(3.157) 

QBA(x)ix)  =  f  k{x,t)-xit)dt, 

Ja 

(3.158) 

QbbUXx)  =  xix)  +  k(x,t)-xit)dt. 

(3.159) 

Given  interpolation  operators 

:  R’>^L^[a,c], 

i^A  :  R^-^L^A), 

0s  :  R'’-I"(B), 

defined  by  (3.14),  we  will  in  this  section  refer  to  them  by  03,  03^, 
The  operator  03  is  related  to  03^ ,  03^  via  the  expressions 

and  03g,  respectively. 

=  ^’3a  ■  C'd  ’ 

(3.160) 

0'3|b  —  V3b  ■  Cu  ) 

(3.161) 

with  Cd,Cu  :  RP  — *  Rp  given  by  (3.21),  (3.22),  respectively.  We 

in  addition  require  the 

zero  operators 


0^,  :  R^(^-^>xL^(Af)^L^(B), 

Oa,  :  R^^L^(B), 

Ob,  :  X  L^(Bf)  -  L^(A), 

Ob,  :  R’^-^L^iA). 

Given  Discretization  D.  we  require  the  interpolation  operators 

i;.-!  :  R”"  X  L^(Bf)  ^  L^la,c], 
02  :  L^iAfjxR^o  ^  L^la,cl 
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given  by  the  formulae 


=  X  0^2  » 

(3.162) 

^1|B 

=  Ob,  X  V’ib  , 

(3.163) 

=  ^2^  X  Oa,  , 

(3.164) 

^^218 

=  0b2  X  V>3g  . 

(3.165) 

Assuming  that  the  operators  Paa^  Pbb  are  nonsingular,  we  define  the  functions 

TfA  ■■  A  —  R, 

7?b  :  B  ~^K, 

via  the  formulae 

VA  =  PAif\A),  (3.166) 

VB  =  PEUf\B)-  (3.167) 

Similarly,  assuming  that  the  operators  Q,QaAiQbb  are  nonsingular,  we  then  define  the 
operators 

XI  :  X  LHBf) L^[a,c]  , 

X2  ■■  L\Aj)xR^'^  -^L^[a,c], 

X3  :  RP-l2[a,c], 

4>u  ■  X  I2(A,) , 

4>2^  :  L\Aj)xR^<^-^^  L^{A)  , 

<h^  ;  RP-XV), 

xL^{Bj)^  L\B)  , 

<hs  ■■  X''(Rc)  X  ^  X2(5)  , 

<hs  ■  HJ>-L\B), 

via  the  formulae 


Xi  =  Q"*(t/’i)  , 

(3.168) 

X2  =  <3"HV’2)  , 

(3.169) 

X3  =  Q~^{^3)  , 

(3.170) 

(3.171) 

=  Q~A\{'i’2y,)  , 

(3.172) 

=  Q'AAi^^A)  ' 

(3.173) 

‘^Ib  -  Qsb(^’Ib)  , 

(3.174) 

—  ^Bb(^2b)  .. 

(3.17.5) 

^  <3bb(^3b)  . 

(3.176) 

1 
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Finally,  given  the  transpose  defined  by  (3.98),  we  will  define  operators 

RP  (?-1)  X  l\Ac)  ^  X  L^(Ac)  , 

0(2 

L^{Aj)  X  RP  (9-i)  ^  RP  lfl-i)  X  X^(Ac)  , 

0^13 

RP  _  RP(</-1)  X  L\Ac)  , 

021 

RP  (9-1)  X  L'^iA^)  ^  L'^iAj)  X  RP  ^’-i)  , 

Q22 

L^{Af)  X  Rp(9-i)  L^{Aj)  X  RP(9-i)  , 

®23 

RP  L^{Aj)  X  RP<’-i)  , 

®31 

Rp.(,-i)  ^  l\Ac)  RP  , 

“32 

L\Aj)  X  RP(9-i)  ^  RP  , 

RP  -  RP  , 

afi 

RP(«-i)  X  L\Bj)  RP-(?-i)  X  L'^iBf)  , 

«f2 

L'^iBc)  X  RP-^«-')  RP(9-i)  X  L'^iBf)  , 

"fs 

RP  RP(?-1)  X  i\Bf)  , 

ofi 

RP(?-i)  X  L\Bf)  ^  L^{Bc)  X  RP(»-i)  , 

a^2 

L^{Bc)  X  RP  (9-i)  ^  L^(Rc)  X  RP  (’-^)  , 

ofs 

RP  -  L\Bc)  X  RP  (?-^)  , 

afi 

Rp.(,-i)  ^  -*  RP  , 

«32 

L^iBc)  X  RP  (<»-J)  -  RP  , 

RP  ^  RP  , 

RP’  X  L'^iBf)  RP’  X  L\Bf)  , 

012 

L^iAj)  X  RP"  -  RP’  X  L\Bf)  , 

Q'lS 

RP  ^  RP’  X  L^{Bj)  , 

021 

RP’  X  L^iBf)  -  L^{Af)  X  RP’  , 

022 

L\Af)  X  RP’  L^(Af)  X  RP’  , 

023 

RP  -  L^iAf)  X  RP’  , 

Q3I 

RP’  X  L'^iBf)  -  RP  , 

O32 

L^(Af)  X  RP’  -  RP  , 

O33 

RP  -  RP  , 
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by  the  formulae 

«2i  =  •  <i>iA  , 

«3i  =  •  <t>ly,  , 

«fi  =  '•PiB  •  ’ 

«fl  =  ^2b  ■  ’ 

=  ‘^Ib  '  ' 

i'l-Xi 
02  -Xl 
03  •  Xl 

^2^ 
^3 

f,B 


Si 

Si 

S3 

V'y,  •  VA 

^'L  • 


ttl2  =  0^ 

•02^  . 

0(3 

=  0u 

•03.4 

II 

'iiS 

0 

■02x  . 

<^23 

=  02^^ 

•03.4 

(^32  =  03^ 

■02^  , 

=  03^ 

•03a 

t>f2  =  0L 

‘02s  * 

«f3 

=  0?B 

•03s 

«f2  =  0Jb 

•  02b  > 

^23 

=  0L 

•  03b 

«f2  =  03^b 

■  02b  ’ 

=  03b 

•  03b 

O12  =  0?" 

•X2  , 

Ol3  = 

=  0r 

X3  , 

022  =  ^2 

•X2  > 

023  = 

=  0^  •: 

X3  > 

O32  =  03 

•X2  , 

O33  = 

=  03  • 

X3  , 

6  RP(9-l) 

X  L\A,)  , 

e  L'^iAj) 

X  RP  <'' 

-1) 

e  R^, 

e  Rp-^’-*) 

X  L\Bf)  , 

e  L\Bc) 

X  RP  ^" 

-1) 

1 

e  R” , 

C  RP  ’  X  L^{B]) 

€  L\Aj) 

X  RP  ’ 

> 

e  R'’ , 

^2  =  02^, 

-^  .4  , 

^3  = 

=  0^- 

'Ha  , 

■riB  . 

= 

=  0Jb  • 

On  = 
0:21  = 

031  = 

and  the  functions 


by  the  expressions 


^i'  = 
= 


(3.177) 


(3.178) 


(3.179) 


(3.180) 


Sy  =  vj  a  , 


Si  -  ^2  •  ' 


^3  =  03  ■  , 
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where  ct  is  the  solution  to  equation  (3.37). 

3.6.2  Recursive  solution  of  the  integral  equation 

We  now  consider  the  originjJ  integral  equation  (3.37) 

Pc  =  f. 

The  main  result  of  this  section  is  the  following  lemma,  which  constructs  the  solution  a 
of  equation  (3.37)  from  the  solutions  tia,t)b  of  equations  (3.166)  and  (3.167). 

Lemma  3.10  If  Discretization  D  is  used,  so  that  operators  Pab,Pba  can  be  approximated 
via  (3.153)-(3.154),  and  if,  in  the  notation  of  Section  3.5.1,  all  six  operators  P,  Paa,  Pbb, 
Q,  Qaa,  Qbb  ore  non-singular,  and  the  operators 

Ai  ;  X  L'^iAc)  ^  x  L^Ac)  , 

A2  ;  L\Bc)  X  RP  («-i)  L^{Bc)  X  Rp  («-i)  , 


given  by 

Ai  =  !Rpii-i)xL^{Ac)  ~  “n  ■  ^AB  ■  0122  • 

(3.181) 

A2  =  ^L^{Bc)xKp(i-^^  ~  0^22  '  ^BA  •  Ctn  •  MaBi 

(3.182) 

are  also  nonsingular,  then  to  accuracy  0(1/ p'‘~^). 

0\A  =  '74+  •  ^AB  •  -Aj  *  •  (ct22  •  ^BA  ■  -  if)  , 

(3.183) 

^\B  =  t/S  +  02b  •  Mb  A  •  Af*  •  (q(*i  ■  MaB  -bf  -  S^)  . 

(3.184) 

Proof.  Using  definitions  (3.39)-(3.42),  the  integral  equation 

Pa  =  f 


can  be  rewritten  in  the  form  (3.43)-f3.44).  The  approximation  (3.59)  and  (3.60)  for  Pab  and 
Pba-  respectively,  can  then  be  used  to  obtain  an  explicit  solution  to  the  coupled  equations 
(3.43)  and  (3.44)  in  terms  of  the  functions  tja,  ^b  defined 

by  (3.166).  (3.167),  (3.171),  (3.172),  (3.173),  (3.174),  (3.175),  (3.176),  respectively.  Indeed, 
applying  the  operator  P^\  to  equation  (3.43)  and  the  operator  Pgg  to  equation  (3.44).  we 
have 

+  Pa\  ■  PABicr^B)  =  PaHAa),  (3-185) 

^BB  ■  PBAi<f\A)  +  ^\B  —  PbB^Ab)-  (3.186) 

Substituting  the  approximations  (3.-59)  and  (3.60)  into  (3.185)  and  (3.186)  yields  the  for¬ 
mulae 

‘^1-4  +  Paa  ■  ■  t’Jg  •  a|jg  = 


(3.187) 
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^BB  ■  •  Mba  ■  ■  <^\A  +  (^\B  =  VB,  (3.188) 

or 

<^\A  +  4>lA-^^AB-i’2B'^\B  =  nA,  (3.189) 

4>2b  •  ^^BA  •  ■  (^\A  +  cr|B  =  VB,  (3.190) 

where  we  have  used  the  definitions  (3.171),  (3.175)  for  4>ij^  and  respectively.  Now, 
multiplying  (3.190)  by  ■  Mab  ■  V’la  subtracting  it  from  (3.189),  we  obtain 

{Il-^(A)-  (t>\A-  ^^ab  ^BA-i>Jj-(J\A  =  (3.191) 

Va  -  4>Ia  •  ^^ab  ■  ■  VB- 

Similarly,  multiplying  (3.189)  by  4>2g  ■  Mba  ■  4a ^  subtracting  it  from  (3.190)  results  in 
the  equation 


( -  <nB  ■  ^^ba  ■  •  Mab  ■  i’Jg)  ■  (^\b  =  (3.192) 

VB- (ha-  ^^ba  ■  4'Ia  ■  ^-4- 

Due  to  (3.177),  (3.178)  and  (3.180),  we  can  rewrite  these  equations  in  the  form 

(IlHA)  -  •  ‘^^.46  ■  {Cl22  ■  ^^BA  •  V-'T^))  ■  <^\A  =  VA  -  <t>\ a  '  AB  '  (>2  »  (3.193) 

ih-^{B)  -  S>2b  '  ^^ba  •  (“u  •  Mab  •  4'2g))  ■  <^\b  =  Vb  -  4>2b  '  ^ba  ■  ■  (3.194) 

By  application  of  Lemma  2.5,  we  obtain 

^\A  —  {Ia  +  '  ^"^^AB  '  {lL2{Br)xRP<’i-'^  ~  ^22  '  ^BA  '  0:u  '  Mab)  (3.195) 

0:22  •  ■  vf^)  -iVA-  4)Ia-  ^^AB  ■  ^2)' 


<^|B  —  (Ib  +  4>2g  ■  Mba  '  (IiiPiv->)xL^{A^)  ~  ’  ^^ab  '  0:22  '  ^^^ba)  (3.196) 

•  ^^Iab  ■  ^’’23)  ■  iVB  -  <>2b  •  -Mba  •  <^1’ )■ 

I  he  equations  (3.183),  (3.184)  are  now  obtained  from  equations  (3.195),  (3.191)  and  equa¬ 
tions  (3.181 ).  (3.182).  □ 

3.6.3  Further  Analytical  Results 

U  p  now  collect  a  number  of  identities  which  are  necessary  for  the  algorithm  to  be  presented 
in  Section  3.8.  First,  we  apply  Lemma  3.10  to  the  particular  cases  /  =  il'\,  f  —  4'2,  f  =  t'3 
to  obtain  analytical  expressions  for  the  functions  \i,  \2  and  X3  defined  in  equations  (3.168) 
and  (3.169). 
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Corollary  3.9  If  Discretization  D  is  used,  so  that  operators  Pab,Pb.a  can  be  approximated 
via  (3.153)-(3.154).  and  if,  in  the  notation  of  Section  3.5.1,  all  six  operators  P,  Paa,  Pbb, 
Q,  Qaa,  Qbb  are  non-singular,  then  to  accuracy  0{\jp^~^), 


Xl|A  =  (<?^A  +  •  AIaB  •  A2  '  -  «f2  •  ^BA  ■  0^13)  X 

(3.197) 

(-dq^  •  AIaB  •  -"fl)  1 

Xl|B  =  i-4>2B  '  ^BA  ■  ■  Cifs)  ^ 

(3.198) 

(d'lfl  +  d>2s  •  ^^ba  ■  Af  *  -  of,  •  Mab  ■  ofi)  > 

X2|,.,  =  (d>2^  +  d>i^  •  Mab  ■  A2 '  •  0:^2  ‘  ^ba  *  ofz)  x 

(3.199) 

( “d>i^  •  Mab  •  A2  ^  •  oifs)  > 

X2|b  =  (-d)2B  •  •  Ar' -of.)  X 

(3.200) 

{03b  +  02b  •  ^^ba  ■  A7'  •  of,  •  Mab  •  “fs)  » 

X3|^  —  ■  Cd  +  ■  ^^AB  ■  <^2  ^  ■  (“^2  ■  ^^BA  '  Otfz  '  Cd  ~  Oi23  '  ,  (3.201) 

X3]fl  =  <hB  -Cu  +  <hB  ■  ^fsA  ■  •  (qh  •  Mab  •  afs  •  Cu  -  «i3  •  (^d)  ,  (3.202) 

where  the  coefficients  af^  and  af^  are  given  by  equations  (3.177)  and  (3.178). 

Proof.  We  only  prove  the  expressions  (3.197),  (3.198);  the  proofs  for  the  remaining  expres¬ 

sions  (3.199),  (3.200)  and  (3.201),  (3.202)  are  nearly  identical. 

Substituting  in  equations  (3.195),  (3.191)  the  functions  d>i^,<f>ig  defined  by  (3.171), 
(3.174)  for  the  functions  t]a,  t]b  defined  by  (3.166),  (3.167),  and  the  matrices  oii\.Oi22  defined 
by  (3.177),  (3.178)  for  the  vectors  defined  by  (3.180),  we  obtain 

Xl|^  =  ■  hlAB  ■  ^2  ^  ■  ^^22  '  ^^BA  '  V'l^)  '  (3.203) 

(d>u  ~  ■  ^^AB  •of2)> 

'Vi|B  =  {Il^{B)  F  (hs  '  ^^BA  ' '  Oiu  '  hlAB  ■  ^23"^ '  (3.204) 

{02b-  ^b'  ^^ba  •  On). 

The  expressions  (3.197),  (3.198)  are  now'  easily  obtained  from  the  equations  (3.203),  (3.204). 

□ 

We  will  also  require  analytical  expressions  for  the  inner  products  ^i,  ^2-  a^nd  63  defined 
in  (3.180)  in  terms  of  the  restricted  inner  products  6^,  bf,  bf,  63  ,  and  bf  defined  in 
(3.180). 


Corollary  3.10  If  Discretization  D  is  used,  so  that  operators  Pab,Pba  can  be  approx¬ 
imated  via  (3.153)-(3.154)<  and  if,  in  the  notation  of  Section  3.5.1,  all  six  operators  P, 
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Paa,  Pbb,  Q,  Qaa,  Qbb  ore  non-singular,  then  to  accuracy 
^1  =  '’/’f  •  O’  =  X  •  <7j^  +  •  ^IB 

/  ^3^  +  ■  Mab  •  AJ'  •  (af2  •  Mba  ■  -  ^f)  \ 


=  ,  (3.205) 

V  +  Qf2  ■  Mba  •  Ar'  ■  (af,  •  Mab  •  ~  6^)  ) 

62  =  =  X0ai)^-<T|.4  +  (0b2  xV’3B)^-t^|B 

/  6^  +  •  Mab  •  Aj  ‘  •  (q®  •  Mba  ■  \ 

=  ,  (3.206) 

V  ^3®  +  0I2  •  Mba  ■  Afi  •  (a^*,  •  Mab  ■  6f  -  6^)  ) 

^3  —  ■  O  =  (i/’3^  •  Cb)^  •  Oj4  -f  (03g  •  Cb)^  •  0|B 

=  C'd  •  ^3  +  Clr  ■  +  P^D  '  “31  ■  Mab  ■  AJ^  •  (0^2  •  Mba  ■  -  <5f ) 

+cj  ■  a§2  ■  Mba  ■  A-i  •  ■  Mab  ■  ^2  "  ^i)  ■  (3-207) 

Proof.  Multiplying  equation  (3.183)  by  and  and  equation  (3.184)  by  and 
V2f.  -  ^’6  obtain 


—  62  +  ct2i  ■  Mab  ■  Aj  ’  •  (o^  ■  Mba ' ~  ^2  )  ■> 
=  +  ai,  ■  Mab  ■  Aj"'  •  (ofj  •  Mba  ■  -  Sf)  , 
=  Sf  +  ofj  •  Mba  •  Af’  •  ■  Mab  ■  , 


'^^3s-^\b  =  ^3  +<^32- BA- -{at, -Mab -Si -Sf)  . 


(3.208) 

(3.209) 

(3.210) 

(3.211) 


Now,  e.xpressions  (3.205),  (3.206)  and  (3.207)  are  easily  obtained  from  (3.208)-(3.211). 

□ 

Special  cases  of  Corollary  3.10  are  obtained  when  /  =  tpi,  f  =  ti>2,  or  /  =  ^3.  While  the 
objects  in  Corollaries  3.11-3.19  are  different  from  those  of  Corollary  3.10  (Corollary  3.10 
is  concerned  with  the  vectors  ^i,52^3»  while  Corollaries  3.11-3.19  are  concerned  with  the 
matrices  qh.  0121  «i3,  Q'2i5  <>22i  c»23t  0311  <*32?  033),  the  proof  for  each  of  CoroDaries  3.11-3.19 
is  nearly  identical  to  that  of  Corollary  3.10. 

Corollary  3.11  //  Discretization  D  is  used,  so  that  operators  PaBiPba  can  be  approx¬ 
imated  via  (3. 153)-(3.154),  ond  if,  in  the  notation  of  Section  3.5.1,  all  six  operators  P, 
P.AA'  Pbb.  Q.  QaA'  Qbb  ore  non-singular,  then  to  accuracy  0{l/p^~^), 


On  =  ti’i  ■  \i  = 


(3.212) 


RP(?-i)  5^  L^(Bj)  —  RP, 

RP  RP(<7-1)  X  L^{Bf), 

Rp.(,-i)  ^  i^^Bf)  -  RP’(9-J)  X  L\Bf), 
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given  by 


II 

0^  +  031  •  M^b  •  ^2  ^  ■ 

a®2  •  Mba  ■  oi3 

II 

—  031  •  Mab  ■  A2  ^  •  0^1 

— 

—0^  •  Mba  '  Ai  ^  •  0(^3 

ofi  +  ofj  •  Mba  •  • 

On  •  M^b  •  On 

Corollary  3.12  If  Discretization  D  is  used,  so  that  operators  Pab,Pba  can  be  approx¬ 
imated  via  (3.153)-(3.154),  and  if,  in  the  notation  of  Section  3.5.1,  all  six  operators  P, 
Faa,  Pbb,  Q,  Qaa,  Qbb  o,re  non-singular,  then  to  accuracy  0{l/p^~^), 


with 


given  by 


012 


=  •  X2  = 


Mpxi  Mpp 
Mjin  Mjip 


L^(Af)  X  RP(<'-i)  RP, 

Mpp  '• 

RP  RP, 

\f 

•‘"nn 

L^(Af)  X  Rp(9-i)  RP(9-i)  X  L^Bj), 

Mfip 

RP  RP(9-i)  X  L^iBf), 

II 

+  «31  •  ^^AB  •  Aj  ^  •  Qf2  •  Mba  ■  o^ 

II 

— 031  •  Mab  ■  A2  ^  •  0^  , 

Mnr.  = 

~Oi2  •  AIba  ■  Aj  ^  •  a^2  » 

Mnp  — 

0^  +  0^  •  Mba  ■  Aj  ^  •  q(\  •  Mab  ■  O23 

(3.213) 


Corollary  3.13  If  Discretization  D  is  used,  so  that  operators  Pab^Pba  can  be  approx¬ 
imated  via  (3.153)-(3.154),  and  if,  in  the  notation  of  Section  3.5.1,  all  six  operators  P, 
Faa-  PbB}  Q-  Qaa-  Qbb  ere  non-singular,  then  to  accuracy  0{l/p^~^), 


with 


given  by 


Mpp 


013  =  ■  X3  = 


Mpp  \ 

Mnp  )  ' 


(3.214) 


Mpp  :  RP-RP, 

A/„p  :  RP  -  RP  (’-')  X  1^(5/), 


o^  •  Cp  +  a3j  •  M ab  ■  ^2  *  ’  ■  ^ttA  ■  0(^3  •  Cp  O23  •  Cp)  , 

ofa  ■  Cu  +  ofj  •  Mba  ■  ■  (®n  ’  ■  0^3  •  Cu  -  0(^3  ■  Cp)  . 
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Corollary  3.14  If  Discretization  D  is  used,  so  that  operators  PaB^Pba  can  be  approx¬ 
imated  via  (3.153)-(3.154),  and  if,  in  the  notation  of  Section  3.5.1,  all  six  operators  P, 
Paa,  Pbb,  Q,  Qaa,  Qbb  are  non-singular,  then  to  accuracy  0(1/?*""^), 


with 


given  by 


021 


=  V’J  •  Xi 


Mfip  M rifi  I 

Mpp  Mpn  )  ’ 


A/„p  ;  L\Af)xW<‘>-^\ 

A/„„  :  X  L\Bf)-^  L\Aj)xR^-^‘‘-''\ 

A/pp  :  RP-RP, 

A/p„  :  RP  (’-‘)  X  1^(5^)  RP, 


A/„p  =  0^3  +  O21  •  MaB  ■  A2  ^  •  0:^2  ■  ^BA  '  C»13  » 
^fnn  —  ~^2\  ’  ^^AB  ‘  ^2  '  ®21  ’ 

A/pp  =  •  A/bx  ■  ■> 

A/pn  =  +  Oc^2  ■  ^^BA  ■  Aj  '  •  O^j  •  MaB  '  C*21  • 


(3.215) 


Corollary  3.15  If  Discretization  D  is  used,  so  that  operators  PaBiPba  can  be  approx¬ 
imated  via  (3.153)-(3.154),  and  if,  in  the  notation  of  Section  3.5.1,  all  six  operators  P, 
Paa,  Pbb,  Q,  Qaa,  Qbb  are  non-singular,  then  to  accuracy  0(l/p*~^), 


with 


given  by 


022 


T 

=  ^*2  ■  \2  = 


■^^nn  np  | 

^fpn  j 


Mnn  ■■  X  -  I2(.4/)  X 

A/„p  ;  RP -- L^(A/)  X  R?4?-'). 

.\/p„  :  A2(.4/)  X  RP'?-0  _  RP, 

.V/pp  ;  RP  RP. 


A/nn  =  ®22  4"  ®21  ’  AB  '  ^2  ®22  '  ^BA  '  ^12  ' 

Alnp  =  ~*^21  ’  AB  '  A2  •  O23  , 

=  ~^32  '  ^^BA  ■  Aj  ■  0^12  ’ 

A/pp  =  O^  +  0®2  ■  A/b>1  •  A,  *  •  O'^i  •  MaB  '  (^23  ■ 


(3.216) 
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Corollary  3.16  If  Discretization  D  is  used,  so  that  operators  Pab,Pba  can  be  approx¬ 
imated  via  (3.153)-(3.154),  and  if,  in  the  notation  of  Section  3.5.1,  all  six  operators  P, 
Paa,  Pbb,  Q,  Qaa,  Qbb  are  non-singular,  then  to  accuracy  0(l/p^“^), 

with 

A/„p  :  RP  L\Af)xW’-^'^~^\ 

Mpp  :  RP^RP, 

given  by 

Mnp  —  023  ■  "I"  at2i  •  MaB  '  ^2  *  ‘  ^^22  ‘  ^BA  '  0^3  •  Cu  —  •  Cjj^  , 

A/pp  =  •  Cfj  +  0^2  •  Mba  ■  ^  •  Mab  '  ■  U-u  ~  0:^3  •  C^j)  . 

Corollary  3.17  If  Discretization  D  is  used,  so  that  operators  PaBiPba  ^an  be  approx¬ 
imated  via  (3.153)-(3.154),  and  if,  in  the  notation  of  Section  3.5.1,  all  six  operators  P, 
Paa,  Pbb,  Q,  Qaa,  Qbb  are  non-singular,  then  to  accuracy  0(l/p'‘~^), 

tt31  =  ^3  -Xl  = 

{Cp  •  +  Cp  •  031  •  Mab  ■  ^2  ^  ‘  ^22  ‘  ^ba  ■  <^13  (3.218) 

—Cjf  •  q®2  •  Mba  '  ^  ^ 

i-Cj)  •  0:31  •  Mab  ■  Aj'  •  ofi  +  ■  of^  + 

Qu  •  <^32  •  ^^ba  •  •  a(\  •  Mab  •  afi)  • 

Corollary  3.18  If  Discretization  D  is  used,  so  that  operators  Pab,Pba  can  be  approx¬ 
imated  via  (3.153)-(3.154),  and  if,  in  the  notation  of  Section  3.5.1,  all  six  operators  P, 
Paa  -  Pbb-  Q-  Qaa-  Qbb  non-singular,  then  to  accuracy 

032  =  rj  ■  X2  = 

(C  Q  ■  Q32  +  Cj)  •  Q31  •  Mab  ■  A2  •  0^2  ‘  ^^b.a  ■  ^12 

—Cjj  ■  0^2  ■  ^^ba  '  A]  •  Q12)  ^ 

( ~^D  •  03J  •  Mab  ■  A2  *  ■  0^3  +  cj  •  + 

C'?  •  «f2  •  •  Ar’  •  ofi  •  Mab  ■  afa)  •  (3.219) 

Corollary  3.19  If  Discretization  D  is  used,  so  that  operators  Pab,Pba  can  be  approx¬ 
imated  via  (3. 153)-(3. 154)-  and  if,  in  the  notation  of  Section  3.5.1,  all  six  operators  P, 
Paa-  Pbb-  Q-  Qaa-  Qbb  are  non-singular,  then  to  accuracy  0(l/p*“'), 

=  vj  ■  X3  =  (3.220) 

C'd  ■  •  Cd  +  Cj)  ■  Q3J  •  Mab  ■  Aj’  •  (Qf2  •  Mba  ■  0113  ■  Co  -  -Cu) 

+Q/  •  •  Cc  +  €[■  ■  0^2  ■  Mba  ■  A;  ‘  •  (q(\  ■  Mab  ■  afa  '  -  0(^3  •  Co)- 


O' 33 
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Finally,  combining  Lemma  3.10  with  the  expressions  (3.197)-(3.202),  we  have 

Corollary  3.20  Suppose  that  Discretization  D  is  used,  so  that  operators  PaBiPba  can  be 
approximated  via  (3.153)-(3.154),  and  suppose  further  that  in  the  notation  of  Section  3.5.1, 
all  six  operators  P,  PaAi  Pbb,  Q>  QAAy  Qbb  are  non-singular.  Suppose  finally  that  the 
function  F  :  [a,c]  — *  R  is  defined  by  the  formula 

F(x)  =  xi(:r) '  (  ^  X2(x) '  (  ^22  )  ^  ^ 

with 

'^Ui'^225-^3  C 

Ai2  €  X 

A21  €  L\Aj)xR^<‘>-^'>. 

Then  on  the  interval  [a, 6],  to  accuracy  0{l/p'^~^), 

F(x)  =  0i^{x)-pi  +  4>2A^)  ■  F2  +  4>3y^{x)  ■  P3  +  VAi^)  ,  (3.222) 

with 

Ml  6 

M2  €  I2(.4/)  X  RP<''-'), 

M3  €  R*', 

defined  by  the  formulae 

Ml  =  MaB  •  -^2  '  ■  i^22  ■  ^^BA  •  (<lu  •  All  +  0^2  '  ^^21  +  “iS  '  C'd  ’  A3  +  ) 

-ofj  •  Ai2  -  0^3  ■  (A22  +  Cy  ■  A3)  -if), 

M2  =  A21  , 

M3  =  Aii+C£)-A3.  (3.223) 

Similarly,  on  the  interi'ol  [6, c],  to  accuracy  0{\lp^~^), 

F{x)  =  Oig{x)  I'l  F  (hgi^) '  ^2  F  ’nBi^)  y  (3.224) 

with 

i/i  €  RP<’-1)  X  1^(5/), 

U2  e  ^  RP  (<7-1)^ 

M3  €  R^, 


defined  by 

the  formulae 

I'l 

=  Ai2  , 

1/2 

■  Mab 

■  (o^  •  A22  +  C»fi  •  Ai2  +  023  ■  ■  A3  +  ^®) 

4  \  ^ 

-Oj2  •  A21  -  Qi3  • 

(All  + 

Cd  A3)-^^  , 

^'3 

=  A22  +  Cy  ■  A3. 

(3.225) 
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Proof.  Restricting  (3.221)  on  the  subintervals  A,B  of  [a,c],  respectively,  we  have 

(3.226) 

Now,  the  expressions  (3.223),  (3.225)  follow  by  combining  (3.226),  (3.227)  with  (3.183), 
(3.184),  (3.197),  (3.198),  (3.199),  (3.200),  and  comparing  the  resulting  expression  with  the 

expressions  (3.222),  (3.224). 

3.7  Description  of  the  Algorithm  for  Singular  Solutions 

In  this  section,  we  present  a  merging  strategy  for  integral  equations  with  end-point  singu¬ 
larities  using  Discretizations  Cl  and  C2.  We  subdivide  the  whole  interval  [a,c]  into  2M 
subintervals,  where  M  is  a  positive  integer.  The  boundary  points  61,62,  are  defined 

by  the  expressions 


61 

=  a. 

(3.228) 

6. 

=  (2<i<^+') 

(3.229) 

6. 

=  (W  +  2<.<2M) 

(3.230) 

^7M+\ 

=  c. 

(3.231) 

We  now  define  subintervals  Bx,  82,-  ■  ■■,  B2M  by  the  formulae 

5.  2 
fi.  2 

=  (6.,6.+i),  (1<1<A/) 

=  (63M+i-i,63M-i-2-i),  (Af  +  1  <  f  <  2A/) 

(3.232) 

(3.233) 

so  that  Bx.B2  and  Bm+x,Bm+2  are  aU  the  same  length,  and  for  all  t  =  2,...,M,  the 
subintervals  5,  and  are  equivalent  in  length,  and  are  twice  the  length  of  S._i  and  v\e 

also  define  subintervals  A2.  A3 . Am.  Am+2,  Am+3,  •  •  • ,  A2M  (for  notational  convenience. 

we  leave  Am+i  undefined)  by  the  formulae 

A 

.42  =  BiU  B2. 

T.  =  A,_,  UP.,  (3  <  i  <  A/) 

M  +  l  =  BmUPm  +  I, 

A,  =  .4._iUP.,  (A/  +  2<i<2A/) 
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3.7.1  Notation 

Generalizing  the  notation  of  Section  3.5,  we  will  denote  by  Pa,^Eb,  restriction  to  the 
interval  A,,  jB,  of  the  integral  operator  P,  respectively,  so  that 


a{x)A 

J 

r^i+i 

/  k{x,t)  ■  cr{t)  dt, 

'b. 

(1  <  2  <  M) 

(3.234) 

<7(Z)  + 

fbiu+i-t 

/  k{x,t)  ■  a{t)dt, 

'bsM-t-l-i 

(M  +  1  <  i  <  2M) 

(3.235) 

<7(l)  + 

f  ^  k{x,t)-(7{t)dt, 

(2<i<  M) 

(3.236) 

PA.i(T)(x)  = 

(r(x)A 

*^1 

/  k{xj)-ait)dt. 

{M  +  2<i<  2M) 

(3.237) 

for  any  <t  e  €  L'^{B,),  respectively.  Similarly,  we  will  denote  by  Qa,-,Qb,  the 

restriction  to  the  interval  >1,,  5,  of  the  integral  operator  Q,  respectively  so  that 


<3b.(x)(3^)  =  + 

k{xj)-x{t)dt,  il<  i<M) 

(3.238) 

QbAxXx)  =  <t(x)  + 

■  x{t)  dt,  (M  +  1  <  i  <  2M) 

(3.239) 

k{x.t)  x{t)dt,  {2<i<M) 

Jbi 

(3.240) 

QaA\){^)  =  + 

/■^zw+i 

/  k{xj)-x{t)dt,  {M  +  2<i<2M) 

•'Gw+i-i 

(3.241) 

for  any  x  €  L'^{Ai)  x  Rp,  x  €  L‘^(B,)  x  RP,  respectively.  For  each  B^  we  define  the  functions 

VB, 

;  5.  -  R  . 

<^1b. 

:  RP  X  £2(5.)- 1^(5.)  , 

:  RP  -  L^{B,)  , 

as  the  solutions  of  the  equations 

PbX'Hb.)  -  Jb,^ 

(3.242) 

QbX'^Xb,'^  = 

(3.243) 

Q B,{d^B,  )  ~  ' 

(3.244) 

Similarly,  we  define  the  functions 

h.4. 

:  A,-R, 

:  RP  X  L^{A,)  -  L^(A,)  , 

:  RP^L^'l.l.), 
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by  the  formulae 

PA,i(rA.)  =  /a,, 

Qa.{X3a,)  = 

provided  the  operators  P.a,,  Pb,iQa,,Q B,  are  nonsingular. 

For  each  i  =  1, . .  .,2M,  we  define  the  operators  and  functions 


(3.245) 

(3.246) 

(3.247) 


by  the  formulae 


«f( 

:  RP 

X  L\B,)  - 

RP  X  Z 

'\B,)  , 

.  RP 

RP  X  L^(Bi)  , 

^31 

:  RP 

T 

X 

R”  , 

:  RP 

^RP  , 

6  RP 

X  L\B,)  , 

63' 

6  R'’ 

afr- 

=  • 

«31* 

= 

“^Ib.  ’  ^33 

=  ^Sb. 

^x'- 

-  • 

VB,  ,  63' 

=  '^k 

•m,  • 

(3.248) 


For  each  i  =  2, . . . ,  M,  M  +  2, . . .,  2A/  we  define  the  operators  and  functions 


-fl' 

:  RP  X 

lHa,)^ 

R”  X  L\Ai), 

:  RP  - 

R"  X  L\A,\ 

®31 

:  RP  X 

L\A.)  - 

R", 

;  RP  - 

R^ 

e  RPx 

L\A.). 

^3' 

€  RP. 

bv  the  formulae 


a 


11 


‘21 


-  xlT 

a  13 

II 

■ 

II 

■  Xlyt.  ! 

“22 

- 

- 

■  ^2.,, 

II 

<5. 

•  fP.4.  , 

II 

e. 

(3.249) 


Finally,  we  define  operators 


L^B-i) 

L\B,) 


L\Bx)  . 
LHB2)  , 
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L^{B\f+2)  — ^  L^{Bm+i)  , 

L^{Bm+i)  — ►  L^{Bm+2)  ■> 

Z2(P.+,)-*R'’x  Z"(^.)  ,  (i  =  2,. 

.,M-  l,M  +  2,. 

.,2M-1), 

y.  L\AC,  ^  .  (i  =  2,. 

..,M-2,M  +  2,. 

.,2M-1), 

RP  X  L\A2m)  -  R'’  X  L'^iAM)  , 

RP  X  L‘^{Am)  Rp  X  L'^{A2m)  , 

which  are  given  by  (3.86),  (3.87),  (3.90),  (3.91). 

3.7.2  Discretization  of  the  Restricted  Integral  Equations 

Choosing  an  integer  p  >  1,  we  construct  the  p  scaled  Chebyshev  nodes  defined  by 
(3.18)  on  each  of  the  intervals  Bi,  i  =  1,2, We  then  discretize  the  three  integral 
equations  (3.242),  (3.243),  (3.244)  via  the  Nystrom  algorithm  based  on  p— point  Chebyshev 
quadrature. 

Remark  3.9  We  use  the  standard  Chebyshev  quadratures  (see  [16],  [18],  [20])  when  the 
particular  operator  being  discretized  is  of  the  form  rn3(P^B),  m3(/(Byi)  (defined  in  Lemma 
3. ,5)  or  of  the  form  m2(P/is),  m2(PB/t)  (defined  in  Lemmas  3.4).  When  the  operator  is  of 
the  form  Paa-,Pbb^  or  of  the  form  mi{PAB)y‘’^i{PBA)  (defined  in  Lemma  3.3),  then  it  is 
discretized  via  the  singular  Chebyshev  quadrature  rules  described  in  Section  3.4.  □ 

The  resulting  approximations  to  the  functions  t}b,,4>1b  ,  (hs  nodes  tjg  wiU  be 

denoted  by 


^B. 

■^flB,)  ’ 

■  •  -  ’ 

^B, 

respectively.  For  each  interval  /1,(1  <  i  <  M  —  1),  we  do  not  construct  approximations 
for  <7.4,.  \i^  ,  \2^  for  the  entire  interval  but  only  for  the  “rightmost”  subinterval  Pi+i. 
Similarly,  for  each  interval  Ai{M  +  1  <  «  <  2M  —  1),  we  construct  approximations  for 
*^4. '\i.A  ■  \2^  on  for  the  “leftmost”  subinterval  R,>i.  The  resulting  approximations  to 
these  functions  at  the  nodes  are  denoted  by 


ttT 

r" 

11 

Remark  3.10  It  is  well-known  that  the  order 
‘^2., I'  to  the  functions  <3;>2.,, 


(i  =  2,. 

..,M,M  +  2,. 

..,2M). 

(f  =  2,. 

..,M,M  +  2,. 

..,2A/), 

(i  =  2,. 

. . ,  M,  Af  +  2, . 

..,2A/). 

f  convergence  of  the  approximations  p,  ;. 
!, ,  is  p.  Since  au  subsequent  steps  in  the 
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construction  of  tin  approximate  solution  3  to  the  integral  equation  (3.37)  are  analytic,  the 
convergence  rate  of  the  full  algorithm  depends  entirely  on  the  parameter  p.  For  example, 
by  using  16  scaled  Chebyshev  points  on  each  subinterval  at  the  finest  level,  one  obtains  a 
sixteenth  order  method. 

The  parameter  p  also  determines  the  order  of  Chebyshev  polynomial  used  to  approximate 
well-separated  regions  of  the  kernel.  For  example,  choosing  p  =  8  results  in  single  precision 
interpolation  accuracy  the  kernels  of  interest,  while  choosing  p  =  16  results  in  double 
precision  accuracy. 

Thus,  p  determines  both  the  order  of  convergence  (for  example,  eighth  order  convergence 
for  p  =  8),  and  the  maximum  accuracy  (for  example,  single  precision  accuracy  for  p  =  8). 
□ 

The  operators  of  the  form  Mjig,  Mb  a,  (^11,0^13,0131  and  functions  of  the  form  ^i,  Aj  all 
have  the  property  that  they  are  composed  of  a  finite  rank  function  (corresponding  to  Cheby¬ 
shev  approximation  for  well- separated  intervals)  and  an  L^  function  defined  for  an  interval 
not  weU-separated.  The  L^  portion  of  each  of  these  operators  is  discretized  at  p  Cheby¬ 
shev  nodes;  we  will  refer  to  the  discretizations  of  these  operators  as  Mab,^BA,  “11,013, 
“ai'^iiAi,  respectively.  We  will  also  refer  to  the  discretizations  033,  63,  A3,  which  are 
equivalent  to  “331^35  A3,  respectively,  since  the  operators  033,  ^3,  A3  are  finite  dimensional. 


3.7.3  Informal  Description  of  the  Algorithm  for  Singular  Solutions 


We  begin  by  directly  solving  the  integral  equations  (3.242),  (3.243),  (3.244)  on  each  of  the 
subinter\'als  Bi.  The  algorithm  then  proceeds  with  an  upward  sweep  for  computing  the  oper¬ 
ators  afi ,  0(3 ,  03  j ,  ,  and  functions  6^' ,  ^3 and  a  downward  sweep  for  computing  Xf'  6 

and  A3 '  €  R'’.  Using  Corollaries  3.3-3.7,  we  first  obtain  data  of  the  form  from 


and  also  obtain  from  Then, 


using  Lemma  3.8  and  CoroOaries  3.2-3.7,  we  obtain  for  each  z  =  3, . . . ,  A/,  iV/  +  3 _ _  2M 

the  functions  ,  the  operators  and  the  functions  6^'  from  the 

functions  t]b,,  ,  the  operators  and  the  functions  ,6f'. 

The  splitting  process  proceeds  in  the  reverse  order.  First,  since  [a,c]  =  A\f  U  A2M,  we 
can  use  Lemma  3.8  to  construct  Xf'^,Xf^'^  (by  inspection,  A^  is  the  expression  to  w^hich 
<^1/4  ‘s  applied  in  equation  (3.116),  Af  is  the  expression  to  which  d>ig  is  applied  in  equation 
(3.117),  and  A3  =  0,  A^  =  0).  Using  (3.146)  from  CoroUary  3.8,  we  can  immediately  obtain 
the  solution  o^b^  using  the  functions  and  the  vectors  A^*^, 

A^^.  Similarly,  we  obtain  the  solution  using  the  functions  aA2M\B2M-  ’ 

A2b,„,„  a-nd  the  vectors  X^^^ .  CoroUary  3.8  also  pro  ides  the  formula  for  the 

2M  1^2  vf  1.  * 

calculation  of  A(*',A3’  given  A^”*^' , Ag'"^'.  Therefore,  we  compute  the  A^',A^'^‘*^'  in  the 
order  i  =  M  -  1.  A/  —  2, . .  .,2,  and  subsequently  use  equation  (3.146)  from  CoroUary  3.8 
to  determine  the  solutions  o\b,,o>\B!^.^,-  tlien  calculate  four  final  functions,  Af'.A®', 

Af^"*^',  e  RP,  from  the  vectors  A(*^,  A3*,  A3  Using  equation  (3.143)  from 

CoroUary  3.8.  we  obtain  the  solutions 
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To  summarize,  the  algorithm  consists  of  three  parts.  First,  a  sufficiently  fine  subdivision 
62, . . 62Af+i  of  the  interval  [a,c]  is  chosen  so  that,  on  each  of  the  intervals  Bi,  the 
functions  >  ^a,ie  ?  Xi a  » X3a  accurately  represented  by  a  low  order 

Chebyshev  expansion  (the  latter  three  functions  are  not  computed  for  both  1  =  1  and  i  = 
A/  + 1).  On  each  of  the  intervals  Bi,  the  equations  (3.242)-(3.244)  are  solved  (approximately) 
by  direct  inversion  of  the  linear  system  arising  from  a  Nystrom  discretization.  Second, 
the  functions  (rA.^g  ,  Xia,|b  ’  tnatrices  and  vectors  S'/  are  computed  for  i  = 

2. . . . ,  M,  A/  +  2, . . . ,  2M.  The  vectors  ,  A3  * ,  are  computed  for  i  =  M,  M  —  1, . .  .,2  and 
also  for  i  =  2M,2M  -  1,...,A/  +  2,  and  finally  the  vectors  Af , Af The 

desired  function  a  is  then  recovered  on  each  subinterval  B,. 

The  following  is  a  more  detailed  description  of  the  numerical  procedure. 

Algorithm  C 

Comment  [Define  the  computational  grid.] 

Create  2A/  subintervals  on  [a,  c]  by  using  the  sequence  of  boundary  points  61,62,..., 

h.v ■  given  by  (3.228)-(3.231).  Create  subintervals  Bi, . . . ,  defined  by  (3.232)-(3.233), 

and  choose  the  number  p  of  Chebyshev  nodes  on  each  Bi.  Determine  the  locations  of  the 

scaled  Chebyshev  nodes  tie,  1  ^2fl, ,  •  •  • , tpa,  on  each  interval  Bi,  and  evcduate  the  functions 

/,  t’l.  il'2  at  these  nodes,  obtaining  /b,,¥’1s  >  •  Compute  the  discretized  operators 

Afe,B,,  AffijB,,  A/bm+jBm+j- A/a.B.+i  .  A/b.+,-4.  (i  =  2, . . . ,  M  -  I, 

Af  +  2 . 2A/  -  1),  A/am/1jm  •  A/ajv-4m  ■  using  one-dimensional  and  two-dimensional 

Chebyshev  transforms,  as  appropriate  (see  Theorem  3.2  and  Lemmas  3. 3-3. 5). 


Step  1. 

Comment  [Construct  the  approximate  solutions  ^b.  ,<^1a  >  *^3 a  on  each  interval  Bi.j 
do  i  =  1,2,. . .,  2M 

(1)  Construct  the  two  p  x  p  matrices  on  B,  obtained  through  a  Nystrom  discretization 
of  the  corresponding  integral  equation. 

(2)  Solve  the  three  p  x  p  linear  systems  on  i?[  by  Gaussian  elimination,  obtaining 
the  values  t/b,  ■  Oib,  .  ^3b, 

end  do 


Step  2. 

Comment  [Construct  the  four  matrices  and  the  two  vectors  on  each 

interval  B, .] 

do  i  =  1,2,. . .,  2M 

Evaluate  the  four  matrices  <5®’ ,  df3’,  Qjj',  d^’  and  the  two  vectors  using  the 

p-point  Chebyshev  quadrature  formula  (see  Remark  3.9). 
end  do 
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jtep  3  (UpwMd  Sweep). 


Comment  [Construct  the  three  functions  ^a.b  iXIa  ^  matrices  a‘^{,  dfs,  031'. 

a^‘  and  the  two  vectors  6^' ,  6^‘  for  i  =  2, . . . ,  M,  M  +  2, ... ,  2M .] 

(1)  Compute  the  three  functions  ,xia  iX3a  >  the  four  matrices  Qn’ , diV , dgi’ , , 

-4-4  ’  -  BBS 

and  the  two  vectors  from  the  data  Mb.Bj, -WbjBi,  «ii' .  “13  >  “3i‘ 

d^‘,  Qf[‘,Qf3^  Qai^a^,  using  the  results  of  Lemma 3.8  and 

Corollaries  3. 2-3. 7.  Similarly,  compute  the  three  functions  >X1aw.  - 

X3^  ,  the  four  matrices  ,  and  the  two  vectors 


from  the  data  Mbs 


+  1  Bm+2  >  ^^M  +  3  ?  ^ 

+  l  +  l  ;r^A/  +  ?  ;r^M+2  -^M+a  C^M  +  l  cBm^2  cBm^^ 

^31  1  “33  ’^11  ’^13  >  *^33  ’^1  » ^3  ’  *^3 

(2)  do  i=  3,4,. .  .,M,M+3,M+4v  .  .,2M 

Compute  the  three  functions  »  Xu  *  X3a  >  matri^'es  af{ ,  af^ ,  Q31 , 0:33 , 

i“»  *1®*  "'>4^4 

and  the  two  vectors  ^[*‘,^3  '  from  the  data  M  j  ,_^b,,  Ms,a,.i  ^b.  .  .  ^3s,  , 

®3i  ”*  ■  ■  ‘^fi '  ®3i‘''^S} '  ^1*”'  ^3  ”* .  ^f'l  '  >  using  the  results  of  Lemma  3.8 

and  Corollaries  3.2 -3.7. 
end  do 


Bm^i  xBm  +  i 


'13 


Step  4  (Downward  Sweep). 

Comment  [Construct  the  two  vectors  AJ*’,  A3’  for  adl  intervals  Ai.  Compute  the  vectors  Af ' ,  A^', 

Af"'*‘,  Af”*‘.] 

(1)  Set  Aj  =  0,  A3 =  0.  Use  the  results  of  Lemma  3.8  to  construct  Xf'^,  Xf^'^  from  the 

data  ttf« ,  ,  6['« ,  . 

(2)  do  I=2M-l,2M-2,  ,2M.M-2,M-3 _ 2 

fjse  Corollary  3.8  to  compute  the  vectors  Af’,  A3  ’  from  the  vectors  lambda'^'*' .  A3 
end  do 

(3)  Use  Corollary  3.8  to  compute  the  vectors  Af'.Aj'  from  the  vectors  Af’.Ag’,  and  tc 
compute  the  vectors  xf'^^'.Xj’^*'  from  the  vectors  A^'*'"*’’ ,  A3 


Step  5. 

Comment  [Compute  the  solution  a  of  equation  (5.37)  at  the  nodes  tis,  ■  ^2b.  ’  ■  ■  ■  1  ^pb, 
interval  Bi.] 

(1)  Determine  the  values  cf  the  the  solution  (T\Bi’<^\Bm+\  using  equation  (3.144). 

(2)  do  1=2,3,.  .  .,M,M+2,M+3 _ 2M 

Determine  the  values  of  the  solution  <r|B,  using  equation  (3.145). 
end  do 
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Remark  3.11  Inspection  of  the  above  algorithm  shows  that  the  amount  of  work  required 
is  of  the  order  0{M-p^).  Step  1  involves  solving  three  pxp  linear  systems  for  each  of  the  2M 
intervals.  Step  2  requires  for  each  interval  the  application  of  two  pxp  matrices  to  two  pxp 
matrices  each  (a  total  of  four  matrix-matrix  multiplies),  and  also  requires  the  application  of 
the  same  two  pxp  matrices  to  a  length  p  vector.  Steps  3-5  require  no  more  than  0{M  -p^) 
operations.  Since  N  =  2M  •  p  is  the  total  number  of  nodes  in  the  discretization  of  the 
interval  [a,c],  we  can  write  the  CPU  time  estimate  in  the  form  0{N  -p^).  □ 

3.8  Description  of  the  Algorithm  for  Smooth  Solutions 

We  turn  now  to  the  construction  of  the  fast  algorithm  for  the  solution  of  the  integral 
equation  (3.37) 

Pc  =  f  , 

using  Discretization  D,  and  based  on  the  apparatus  developed  in  Section  3.6.  The  main 
tool  at  our  disposal  is  the  ability  to  merge  the  solutions  of  restricted  versions  of  the  integral 
equation  in  adjacent  subintervals  (Lemma  3.10).  As  this  suggests  a  recursive  procedure,  we 
begin  by  subdividing  the  whole  interval  [a,c],  on  which  the  solution  to  (3.37)  is  sought,  into 
a  large  number  of  subintervals.  For  the  sake  of  simplicity,  we  assume  that  m  is  a  positive 
integer  and  that  M  =  2”*  is  the  number  of  subintervals  created.  The  boundary  points  of 
the  subintervals  are  then  defined  by  a  strictly  increasing  sequence  of  numbers 

bi,b2,...,bM,bM+u  (3.250) 

with  6i  =  a  and  6Af+i  =  c.  For  each  i  =  1,...,M,  we  define  the  interval  BJ"  via  the 
expression 

Br  =  [bi,bi+r],  (3.251) 

and  create  a  hierarchy  of  intervals  by  recursively  merging  adjacent  pairs.  That  is,  for 
each  j  =  m  -  1, . . . ,  1, 0,  and  i  =  1, . . . ,  M,  we  define 

=  (3-252) 

We  will  refer  to  each  fixed  /  as  a  level.  We  will  also  refer  to  the  two  intervals  and 

as  children  and  to  the  larger  interval  B\  as  a  parent. 

It  is  obvious  that 

»^i-K-2”>-']»  (3.253) 

and  that  for  each  level  /, 

2' 

[a,c]=UB'.  (3.254) 

i=i 

Furthermore,  since  in  this  section  we  are  using  discretization  D,  so  that  the  number  of 
points  in  a  subinterval  is  proportional  to  the  size  of  the  interval),  the  dimension  of  the  finite 
rank  portions  of  Mab-,Mba  given  by  (3.61),  (3.61),  for  an  interval  B\  is  p  •  (m  -  /). 
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3.8.1  Notation 


Generalizing  the  notation  of  Section  3.6,  we  will  denote  by  Pij  the  restriction 
Bi  of  the  integral  operator  P,  so  that 


to  the  interval 


<t(x)+  /’■*■**  k{x,t)  •  <T{t)dt 

(3.255) 

for  any  a  6  L^{B[).  Similarly,  we 

®l+(i-l)  2'"-' 

will  denote  by  Qij  the  restriction  to  the  interval  B\  of 

the  integral  operator  Q,  so  that 

x(x)+/*'^‘*  k{x,t)-xit)di 

(3.256) 

for  any  x  €  L^{B[)  x  Rp  I'^'O.  For  each  B[  we  will  define  the  functions 

r)i.l  : 

B\  ->  R, 

(3.257) 

Rp.(m-0  ^  £^(R-), 

(3.258) 

: 

L\B\)  X  RP-("‘-')  ^  L‘^{B\), 

(3.259) 

W>^L\B[), 

(3.260) 

as  the  solutions  of  the  equations 

PiAVij)  =  f\Bh 

(3.261) 

(3.262) 

(3.263) 

—  ^3|g(  ? 

(3.264) 

provided  the  operators  Pi^uQi^i  are  nonsingular. 

For  each  /  =  0, 1, . .  .,m,  and  i  =  1,2, . .  .,2^  we  define  the  operators 


a 


a 


a 


Q 


Q 


a 


a 


a 


Q 


11 

1.1 
12 

«./ 

13 

1.1 
21 

»,/ 

22 

i,l 

23 

A, 

31 
A, 

32 
A, 

33 


Rp  (m-/)  ^  RP  (m-0  L^(bI), 

L\bI)  X  RP  (”*-0  X  l\bI), 

RP  _  RP  (’"-0  X  l\bI), 

RP  (m-o  ^  -*  L'^{B[)  X  RP  ("-'), 

L^{Bl)  X  RP  <”-')  ^  L\bI)  X  RP-^'"-'), 
L^(B‘i)  X  RP  -  RP  ("*-'), 

Rp  (m-0  ^  ^  RP^ 

L\Bi)  X  RP  ^”*-')  ^  RP, 

RP  ^  RP, 
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by  the  formulae 


«i’l  = 

Qj2 

=  <. 

■4>2i,i  , 

iyl 

“13  = 

Oji  = 

022 

= 

‘^i,l  > 

“23  = 

(3.265) 

«31  = 

^.,1  > 

i,l 

«32 

II 

“33  = 

and  the  functions 

€ 

RP(»n-' 

')  X  L\B^^, 

4'‘ 

€ 

L^Bi) 

X  RP  ("‘-0, 

63 

€ 

R”, 

bv  ♦  he  formulae 

=  ^>1 

.1  '  Vi,l  , 

II 

•Vi,l  , 

^'3  =  ’ 

Pl,  • 

Vi,i  • 

(3.266) 

Finally  for  each  1  = 

=  0,1, 

. . . ,  m  ■ 

-1, 

and  t  = 

1,2,.., 

,2',  we 

define  the  operators 

"2i  — l"2i 

Mg.+ig.+i  :  R'’(’”-')xZ;2(B'tij)^i2(B'f)xR^(”-'), 
defined  by  (3.61)  and  (3.62). 


3.8.2  Discretization  of  the  Restricted  Integral  Equations 

Choosing  an  integer  p>  1,  we  construct  the  p  scaled  Chebyshev  nodes  defined  by  (3.18) 

i 

on  each  of  the  intervals  t  =  1,2,...,M.  We  then  discretize  the  two  integral  equations 
(3.261),  (3.264)  via  the  Nystrom  algorithm  based  on  p-point  Chebyshev  quadrature  (see 
Remark  3.9).  The  resulting  approximations  to  the  functions  p,,/,  It® 

nodes  rf  will  be  denoted  by 

Vi, I 

^.,1 
^,.1 

respectively. 

The  operators  of  the  form  Mab,  Mba,0‘u,oi2,  OiSjOzi)  022,023,  031,032  and  functions 
of  the  form  61 , 62,  -^1 ,  -^2  all  have  the  property  that  they  are  composed  of  a  finite  rank  function 


~  {Vi,hVi,h  •  •  •  iVi,!)  ■> 

=  (*^2,.,>^2,,i’- •  ’ 
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(corresponding  to  Chebyshev  approximation  for  well-separated  intervals)  and  an  L^  function 
defined  for  an  interval  not  weU-separated.  The  portion  of  each  of  these  operators  is 
discretized  at  p  Chebyshev  nodes;  we  will  refer  to  the  discretizations  of  these  operators  as 
Mab  ,  Mb  A,  an ,  "12,  ai3,  021 , 022, 023, 031, 032,  “33,  ^1 ,  ^2,  ,  A2,  respectively.  We  will  also 

refer  to  the  discretizations  033,  ^3,  A3,  which  are  equivalent  to  033,^3,  A3,  respectively,  since 
the  operators  033,  S3,  A3  are  finite  dimensional. 


3.8.3  Informal  Description  of  the  Algorithm  for  Smooth  Solutions 


We  begin  by  directly  solving  the  two  integral  equations  (3.261),  (3.264)  on  each  subinterval 
at  the  finest  level,  as  discussed  in  the  preceding  section.  Theorem  (3.20)  then  shows  that 
<7  restricted  to  can  be  expressed  as  a  linear  combination  of  the  four  solutions  <^ii 

Thus,  it  remains  only  to  determine  the  vectors  Aj”*,  Aj"*  ^  g 

RP  for  each  of  the  M  subintervals  B^.  Fortunately,  this  can  be  done  recursively.  To  see 
this,  suppose  that,  at  some  coarse  level  /  <  m  —  1,  we  are  given  the  vectors  A’/,  Aj^  X'3 
for  the  subinterval  5,-.  Then  Corollary  3.20  provides  formulae  for  the  calculation  of  the 
corresponding  vectors  A^ Af A^‘’'+^  €  € 

RP  for  the  two  child  intervals  B^}^  and  respectively.  On  the  coarsest  level,  we 


observe  that  A°’^  =  0,  Aj’^  =  0,  A®’^  =  0,  i.e.  the  solution  of  equation  (3.261)  on  the  whole 
interval  [a,  c]  is  simply  a. 

However,  the  formulae  (3.223)  and  (3.225)  of  Corollary  3.6  contain  the  eighteen  matrices 
o*t~^’*^^,  (1  <  <  3)  and  the  six  vectors  ^^*’^''■^(1  <  s  <  3).  These  quan¬ 

tities  are  also  computed  recursively  but  in  the  opposite  direction,  namely,  from  the  finest 
level  to  the  coarsest.  They  are  certainly  available  at  level  m  directly  from  the  definitions 
(3.265)-(3.266).  For  the  interval  Bl  at  any  coarser  level  /  <  m  —  1,  Corollary  3.10  and 


Corollaries  3.11-3.19  describe  how  the  nine  matrices  d’’/(l  <  5,t  <  3)  and  the  three  vectors 
^i’^(l  <  5  <  3)  are  obtained  from  the  eighteen  matrices  q:,<,  and  six  vectors  S3  of  the  two 
child  intervals. 

To  summarize,  the  algorithm  consists  of  three  parts.  First,  a  sufficiently  fine  subdivi¬ 
sion  61 , 62,  •  •  • ,  ^M+i  of  the  interval  [a,  c]  is  chosen  so  that,  on  each  of  the  intervals  ,  the 
functions  i/i,  accurately  represented  by  a  low  order  Chebyshev 

expansion.  On  each  of  the  intervals  5,-,^,  the  equations  (3.261)-(3.264)  are  solved  (ap¬ 
proximately)  by  direct  inversion  of  the  linear  system  arising  from  a  Nystrom  discretization. 


Second,  the  matrices  a’’/  and  three  vectors  S'/  are  computed  in  an  upward  sweep,  beginning 
at  the  finest  level  m.  Finally,  the  vectors  Aj^,  x/,  X/  are  computed  in  a  downward  sweep, 
beginning  at  the  coarsest  level.  The  desired  function  <7  is  then  recovered  on  each  subinterval 
from  equations  (3.221)-(3.222). 

The  following  is  a  more  detailed  description  of  the  numerical  procedure. 
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Algorithm  D 

Comment  [Define  the  computational  grid.] 

Create  M  ■=  2'"  subintervals  on  [a,  c]  by  choosing  a  sequence  of  boundary  points  h\,h2, . . . , 

1  with  bi  =  a  and  6*/+!  =  c.  Choose  the  number  p  of  Chebyshev  nodes  on  each 
interval  BJ"  =  [6i,  6,+i]  for  i  =  1, . . . ,  M.  Determine  the  locations  of  the  seeded  Chebyshev 
nodes  , . . . ,  rf  on  each  interval  B[" ,  and  evaluate  the  functions  /,  V'l ,  ti>2 ,  V’s  at  these 
nodes,  obtaining /i_^, For  each  /  =  0,l,...,m-  1,  and  i=  1,2,. ..,2', 
compute  the  discretized  operators  ,  using  one-dimensional  and 

two-dimensional  Chebyshev  transforms,  as  appropriate  (see  Theorem  3.2  and  Lemmas  3.3-3.5). 


Step  1. 

Comment  [Construct  the  approximate  solutions  ^2.,„,^3.,„  on  each  interval  B[".] 

do  i  =  1,2,. . .,  M 

(1)  Construct  the  two  pxp  matrices  on  B|  obtained  through  a  Nystrom 
discretization  of  the  corresponding  integral  equation. 

(2)  Solve  the  four  pxp  linear  systems  on  B[  by  Gaussian  elimination, 
obtaining  the  values 

end  do 


Step  2. 

Comment  [Construct  the  nine  matrices  5*’/"  and  six  vectors  5)  '”  on  each  interval  B["  at  the  finest 
level.] 

do  i  =  1,2,. . .,  M 

Evaluate  the  nine  matrices  dt|'/"  and  three  vectors  5^’"*  using  p— point 
Chebyshev  quadrature  formula  (see  Remark  3.9). 
end  do 


Step  3  (Upward  Sweep). 


Comment  [Construct  the  matrices  d';/  and  vectors  i*  '  for  all  intervals  at  all  coarser  levek  i  = 

m-  l,m- 2,..,,0.] 

do  1=  m-1,  0,  -1 
do  i=l,  a’ 


Compute  the  nine  matrices  q,’,  and  three  vectors  from  the  corresponding  data  in 
the  two  child  intervals  (A/  *■'  -2t-i./+i  .2t  i4.i  ii.,  lo.  i,,. 


B'+'  b'"*’' 


.M 


using  the  results  of  Corollaries  3.4  and  3^5 
end  do 
end  do 


b1+*b^+_*,  ’ 


-2.-l,/+l  -2<,l+l 
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Step  4  (Downward  Sweep). 

Comment  [Construct  the  three  vectors  AJ’”*  for  all  intervals  at  the  finest  level.] 

Set  A°'^  =  0,  X°/  =  0,  Ag’^  =  0. 

do  1=0^-1 
do  i=l,  2‘ 

Use  Corollary  3.20  to  compute  the  vectors  Aj‘’'^^(l  <  s  <  3), 

for  the  child  intervals  and  from  the  vectors  Ai’*(l  <  s  <  3)  of  the  parent 

interval  B[. 
end  do 
end  do 


Step  5. 


Comment  [Compute  the  solution  a  of  equation  (3.37)  at  the  nodes  . .  .,tpB,n  for  each 

interval  BJ"  at  the  finest  level.] 

do  i=l,  M 
doj=l,p 

Determine  the  values  of  the  solution  <r  of  equation  (3.37)  at  the  node  via 
formulae  (3.221)-(3.222). 
end  do 
end  do 


Remark  3.12  Inspection  of  the  above  algorithm  shows  that  the  amount  of  work  required 
is  of  the  order  0{M  •  p^).  Step  1  involves  solving  four  p  X  p  linear  systems  for  each  of  the 
M  intervals.  Step  2  requires  for  each  interval  the  application  of  three  p  x  p  matrices  to 
three  p  x  p  matrices  each  (a  total  of  nine  matrix-matrix  multiplies),  and  also  requires  the 
application  of  the  same  three  px  p  matrices  to  a  length  p  vector.  In  step  3,  the  asymptotic 
cost  for  each  B,-  is  bounded  by  the  cost  of  multiplying  the  largest  matrices  from  the  previous 
level  (1  <  <  3),  each  of  which  have  dimensions  p-(m-/)  xp-(m-/)). 

The  asymptotic  total  cost  is  given  by  the  series 

m— 1 

^  p3  •  2'-*  ‘  •  (m  -  /)  =  4  •  p^  •  (2"*  -  1)  -  2  •  p3  .  m.  (3.267) 

1=0 


Steps  4-5  require  no  more  than  0{M  -p^)  operations.  Since  N  =  M  -  pis  the  total  number 
of  nodes  in  the  discretization  of  the  interval  [a,c],  we  can  write  the  CPU  time  estimate  in 
the  form  0(N  -p^). 
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3.9  Numerical  Results 

FORTRAN  programs  have  been  written  implementing  the  algorithms  described  in  this 
chapter.  In  this  section,  we  discuss  several  details  of  our  implementation,  and  demonstrate 
the  performance  of  the  scheme  with  numerical  examples. 

The  following  technical  details  of  our  implementation  appear  to  be  worth  mentioning. 

1.  Algorithm  C  depends  for  its  stability  on  (3.242)-(3.247)  having  unique  solutions  for 
all  subintervaJs  Ai  and  B,,  while  Algorithm  D  depends  on  the  equations  (3.261)-(3.264) 
having  unique  solutions  for  all  subintervals  Bj  (/  =  0,1,..., M,  i  =  1,...,2^).  It  is  easy 
to  construct  examples  for  which  these  conditions  are  violated,  even  though  equation  (3.37) 
has  a  unique  solution.  In  such  cases,  a  different  subdivision  of  the  interval  [a,c]  can  be 
attempted,  such  that  none  of  the  subintervals  of  the  new  subdivision  coincides  with  an 
interval  of  the  ori^nal  one.  This  procedure  can  be  viewed  as  a  form  of  pivoting,  and  it 
is  easy  to  show  that  it  is  always  possible  to  make  it  work.  However,  since  the  numerical 
ranks  of  the  discretized  operators  of  the  form  Pj^b  ot  Pba  sensitive  to  the  sizes  of  the 
subintervals  A,  B,  an  implementation  of  this  pivoting  scheme  would  be  somewhat  involved. 
It  has  not  been  implemented  at  this  point,  and  we  have  not  so  far  encountered  a  need  for 
it. 

2.  We  have,  however,  implemented  a  crude  scheme  for  detecting  high  condition  numbers  in 
the  algorithms.  For  Algorithm  C,  these  can  occur  in  the  solution  of  the  linear  systems  on 
each  of  the  subintervals  B,  (Step  1),  while  computing  inverses  of  the  matrices  A:  defined  by 
(3. ’30)  used  when  merging  solutions  (Step  3),  and  while  computing  the  inverse  Aj  defined 
by  (3.115)  (Step  4.1).  For  Algorithm  D,  these  can  occur  in  the  solution  of  the  linear  systems 
on  each  of  the  finest  level  subintervals  (Step  1),  and  while  computing  inverses  of  the  matrices 
Ai,  A2  defined  by  (3.181),  (3.182)  (Step  3).  In  all  cases,  the  condition  number  of  the  system 
being  solved  is  estimated  in  the  process  of  solution  (we  use  a  standard  LINPACK  routine), 
and  the  largest  of  these  is  returned  to  the  user.  When  an  extremely  large  condition  number 
is  detected  by  the  LINPACK  routine,  the  resulting  solution  of  the  original  integral  equation 
should  be  viewed  as  suspect.  It  is  easy  to  show  that  when  the  differential  operator  is  positive 
definite,  this  cannot  happen.  A  more  complete  treatment  of  this  subject  requires  further 
study. 

The  algorithms  of  this  chapter  have  been  applied  to  a  variety  of  problems.  Five  experi¬ 
ments  are  described  below,  and  their  results  are  summarized  in  Tables  3.2-3.14. 

In  each  of  these  tables,  the  first  column  contains  the  total  number  N  of  nodes  in  the 
discretization  of  the  interval  [a,  c].  The  second  column  contains  the  relative  L^  error  of 
the  numerical  solution  as  compared  with  the  analytically  obtained  one  at  5000  equispaced 
points  within  the  interval  [a,c],  where  Chebyshev  interpolation  has  been  used  to  evaluate 
the  numerical  solution  at  each  of  the  5000  points.  The  third  column  contains  the  maximum 
absolute  error  obtained  at  any  of  the  5000  points.  The  fourth  column  contains  the  CPU 
time  required  to  solve  the  problem,  excluding  the  time  used  to  evaluate  the  solution  at 
5000  equispaced  points,  where  in  all  cases  the  times  are  given  for  a  Sun  SPARCstation  2 
computer. 
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Table  3.1:  LU  Factorization  timings 


n 

t  (sec.) 

20 

0.100  xlO-^ 

40 

0.300  xlO"^ 

80 

0.180  xl0° 

160 

0.127  xlOi 

320 

0.948  xlO^ 

640 

0.723  xlO^ 

1280 

0.578  xlO^ 

2560 

0.462  xlO"  (est) 

Table  3.2:  Numerical  results  for  Example  3.1,  p  =  5. 


n 

E°-(ct) 

t  (sec.) 

20 

0.184  xl0“ 

0.952  xl0“ 

0.500  xlO'^ 

40 

0.112  xl0° 

0.687  xlO® 

0.150  xl0° 

80 

0.260  xlO-* 

0.152  xlO® 

0.260  xl0° 

160 

0.292  xlO-^ 

0.125  xlO-2 

0.510  xl0° 

320 

0.241  xlO-3 

0.471  xlO-3 

0.106  xlO^ 

Remark  3.13  For  comparison,  Table  3.1  shows  times  for  solving  linear  systems  of  com¬ 
parable  sizes  using  the  LINPACK  LU  factorization  routines  on  a  Sun  SPARCstation  2 
computer.  □ 


Example  3.1  The  following  problem  was  presented  in  [24]  for  the  purpose  of  modeling  the 
intrinsic  viscosities  of  flexible  macromolecules.  The  equation  to  be  solved  is  given  by  the 
formulae 

The  closed  form  solution  of  this  problem  (due  to  [7])  is  given  by  the  formula 


<7(x) 


TT  •  \/2  ■  (1  -  X^)* 


(3.269) 


We  solve  (3.268)  using  Algorithm  C  with  the  number  of  Chebyshev  nodes  p  =  5, 10,20.  The 
results  of  this  experiment  are  presented  in  Tables  3.2-3.4. 


Example  3.2  The  following  problem  was  presented  in  [12]  and  [31]  (in  a  slightly  different 
form)  for  the  purpose  of  modeling  an  antiplane  elasticity  problem  of  a  crack  terminating 
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Table  3.3:  Numerical  results  for  Example  3.1,  p  =  10. 


n 

£V) 

E°°ia) 

i  {sec.) 

40 

0.112  xlO" 

0.686  xlO* 

0.230  xl0° 

80 

0.540  xl0-‘ 

0.335  xl0° 

0.550  xl0° 

160 

0.207  xlO-^ 

0.166  xl0° 

0.123  xlO^ 

320 

0.346  xlO-^ 

0.240  xlO-® 

0.252  xlO^ 

640 

0.789  xlO’* 

0.213  xlO"^ 

0.508  xlO^ 

Table  3.4:  Numerical  results  for  Example  3.1,  p  =  20. 


n 

E°°{a) 

i  (sec.) 

80 

0.541  xl0-‘ 

0.335  xlO" 

0.138  xlO^ 

160 

0.263  xlO-' 

0.149  xl0“ 

0.343  xlO^ 

320 

0.988  xlO-2 

0.720  xlO"^ 

0.738  xlO^ 

640 

0.345  xl0-» 

0.268  xlO-^ 

0.154  xlO^ 

1280 

0.138  xl0-^° 

0.116  xlO-^ 

0.319  xlO^ 

perpendicularly  at  a  bimaterial  interface  (31).  The  equation  to  be  solved  is  given  by  the 
formulae  ^ 

/  I — ^  +  — 1  =  4a:  —  2\/i  +  x^.  (3.270) 

Jo  t  x  +  ifj 

The  closed  form  solution  of  this  problem  (due  in  part  to  [26])  is  given  by  the  formula 

<r(i)  =  — \/i  -  (3.271) 

TT 

We  solve  (3.270)  using  Algorithm  C  with  the  number  of  Chebyshev  nodes  p  =  5, 10, 20.  The 
results  of  this  experiment  are  presented  in  Tables  3.5-3.7. 


Table  3.5:  Numerical  results  for  Example  3.2,  p  =  5. 


’Ll 

E^(a) 

t  (sec.) 

mi 

0.482  xlO-* 

0.102  xlO*' 

Bi 

0.148  xlO-* 

0.482  xlO-^ 

0.120  xl0° 

80 

0.100  xlO-2 

0.619  xlO-2 

0.230  xlO® 

160 

0.288  xlO"^ 

0.161  xl0“^ 

0.460  xlO® 

320 

0.155  xlO-2 

0.115  xlO-^ 

0.102  xlO^ 
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Table  3.6:  Numerical  results  for  Example  3.2,  p  =  10. 


n 

E“(tT) 

i  (sec.) 

40 

0.129  xlO"^ 

0.469  xlO-^ 

0.240  xl0“ 

80 

0.390  xlO'2 

0.201  xlO-^ 

0.550  xl0° 

160 

0.277  xlO"^ 

0.190  xlO-2 

0.119  xlO^ 

320 

0.550  xlO'® 

0.244  xlO-5 

0.237  xlO^ 

640 

0.181  xlO'® 

0.124  xlO-® 

0.489  xlO^ 

Table  3.7:  Numerical  results  for  Example  3.2,  p  =  20. 


n 

E~(<t) 

t  (sec.) 

80 

0.350  xlO"'" 

0.201  xlO"^ 

0.131  xlO^ 

160 

0.102  xlO"^ 

0.680  xlO-2 

0.336  xlO^ 

320 

0.699  xlO"* 

0.464  xlO-3 

0.728  xlO^ 

640 

0.144  xlO'® 

0.862  xlO-® 

0.151  xlO^ 

1280 

0.854  xlO"*® 

0.564  xlO'® 

0.311  xlO^ 

Example  3.3  This  example  resembles  an  experiment  presented  in  [5].  The  equation  to  be 
solved  solved  is  given  by  the  formulae 

A-ct(i)+  /  [(1  +  sin(25x))log(|x  -  t|)  +  cos(25a:t)]  •  =  (3.272) 

Jo 

A  •  sin(mx)  +  gi(x)  •  (1  +  sin(25x))  +  92(2:), 


with  m  =  250,  A  €  {0, 1},  and  71,92  given  by  the  expressions 

7i(x)  =  —  •  [log(x)  —  cos(m)log(l  -  x)  —  cos(mx) 

m 

[Ci(mx)  —  Ci(m(l  —  x))]  —  sin(mx)[Si(mx)  +  Si(m(l  -  x))]], 
_  cos(25x  +  m)(25x  —  m)  —  cos(25x  -  m)(25x  +  m)  +  2m 
^  — 2(25x  +  m)(25x  —  m) 


(3.273) 

(3.274) 


where  Ci  and  Si  are  the  cosine  integral  and  sine  integral,  respectively.  The  solution  to  this 
equation  is  given  by 

<t(x)  =  sin(mx).  (3.275) 


We  first  solve  (3.272)  setting  A  =  1  (a  second  kind  integral  equation),  applying  Algorithm 
D  to  this  equation  with  the  number  of  Chebyshev  nodes  p  =  5, 10,20.  The  results  of  this 
experiment  are  presented  in  Tables  3.8-3.10.  We  then  solve  (3.272)  with  A  =  0  (a  first  kind 
integral  equation),  applying  Algorithm  D  to  this  equation  with  the  number  of  Chebyshev 
nodes  p  =  10,20  (due  to  the  high  condition  number,  p  =  5  yields  no  accuracy  for  this 
problem).  The  results  of  this  second  experiment  are  presented  in  Tables  3.11-3.12. 
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Table  3.8:  Numerical  results  for  Example  3.3,  A  =  l,p  =  5. 


n 

E°°(a) 

t  (sec.) 

20 

0.728  xl0‘ 

0.142  xlO'-^ 

0.600  xlO-i 

40 

0.271  xlO^ 

0.589  xlO’ 

0.310  xl0° 

80 

0.125  xlO^ 

0.218  xlO^ 

0.920  xl0° 

160 

0.201  xl0° 

0.345  xl0° 

0.228  xlO^ 

320 

0.945  xlO-2 

0.163  xlO-^ 

0.532  xlO^ 

640 

0.121  xlO-2 

0.234  xlO-2 

0.122  xl02 

1280 

0.113  xlO-2 

0.190  xlO-2 

0.259  xlO^ 

2560 

0.112  xlO-2 

0.188  xlO-2 

0.553  xl02 

Table  3.9:  Numerical  results  for  Example  3.3,  A  =  l,p=  10. 


n 

E°^ia) 

i  (sec.) 

40 

0.280  xl0‘ 

0.641  xlO^ 

0.390  xlO" 

80 

0.127  xlO^ 

0.273  xlO^ 

0.159  xlO^ 

160 

0.947  xl0-‘ 

0.138  xlO° 

0.495  xlO^ 

320 

0.231  xlO-3 

0.333  xlO-3 

0.135  xl02 

640 

0.297  xlO"® 

0.539  xlO-® 

0.329  xlO^ 

1280 

0.885  xlO-’’ 

0.146  xlO-® 

0.751  xlO^ 

2560 

0.885  xlO-^ 

0.146  xlO"® 

0.164  xlO^ 

Table  3.10:  Numerical  results  for  Example  3.3,  A  =  l,p  =  20. 


n 

E°^ia) 

t  {sec.) 

80 

0.123  xlO‘~ 

0.217  xlO* 

0.236  xlO* 

160 

0.217  xlO-‘ 

0.354  xlO-‘ 

0.105  xlO^ 

320 

0.2U  xlO"® 

0.288  xlO"® 

0.344  xlO^ 

640 

0.321  X  10-^2 

0.467  X  10-^2 

0.943  xlO^ 

1280 

0.107  xl0-*3 

0.340  X  10-^3 

0.237  xlO^ 

2560 

0.112  xlO“^^ 

0.310  xlO"^^ 

0.558  xlO^ 
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Table  3.11:  Numerical  results  for  Example  3.3,  A  =  0,p  =  10. 


n 

E^{a) 

t  (sec.) 

40 

0.102  xlO^ 

0.165  xlO^ 

0.420  xlO° 

80 

0.927  xl0° 

0.146  xlO^ 

0.171  xlO^ 

160 

0.425  xl0° 

0.640  xlO^ 

0.537  xlO^ 

320 

0.740  xlO-^ 

0.135  xlO"* 

0.143  xlO^ 

640 

0.313  xlO'^ 

0.791  xlO'2 

0.351  xlO^ 

1280 

0.861  xlO"^ 

0.200  xlO-^ 

0.801  xlO^ 

2560 

0.202  xlO-2 

0.414  xlO"^ 

0.174  xlO^ 

Table  3.12:  Numerical  results  for  Example  3.3,  A  =  0,p  =  20. 


n 

E^{ct) 

t  (sec.) 

80 

0.952  xlO'^ 

0.160x10^ 

0.242  xlO* 

160 

0.207  xlO-' 

0.417  xlO'^ 

0.108  xlO^ 

320 

0.150  xlO-® 

0.423  xlO*"* 

0.355  xlO^ 

640 

0.345  xlO-^ 

0.699  xlO-® 

0.978  xlO^ 

1280 

0.200  xlO-’’ 

0.436  xlO-® 

0.245  xlO^ 

2560 

0.330  xlO'® 

0.958  xlO"^ 

0.575  xlO^ 
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The  following  observations  can  be  made  from  Tables  3.1-3.12,  and  are  corroborated  by 
our  more  extensive  experiments. 

1.  For  Algorithm  D,  the  practical  convergence  rate  of  the  method  is  consistent  with  the 
theoretical  one. 

2.  While  we  have  not  analyzed  the  theoretical  convergence  rate  for  Algorithm  C,  in  practice 
the  convergence  is  determined  by  the  number  of  subintervals  used,  as  opposed  to  the  total 
number  of  points  in  the  discretization.  While  double  precision  arithmetic  permits  a  maxi¬ 
mum  of  104  subintervals,  our  experiments  indicate  that  maximum  accuracy  is  achieved  by 
using  approximately  64  subintervals. 

3.  Example  3.3  with  A  =  0  is  an  extremely  ill-conditioned  problem,  which  substantially 
reduces  the  accuracy  of  the  computed  results  compared  to  the  relatively  well-conditioned 
Example  3.3  with  A  =  1.  Because  Algorithm  D  is  a  direct  method,  the  timings  for  A  =  0 
and  for  A  =  1  are  the  same,  for  equivalent  number  of  points  n  and  order  of  method  p. 

4.  For  both  Algorithm  C  and  Algorithm  D,  most  of  the  computational  effort  is  devoted 
to  merging  the  a  matrices  (Step  3  in  both  Algorithm  C  and  D).  However,  the  size  of  the 
matrices  is  fixed  for  Algorithm  C,  while  for  Algorithm  D  the  matrices  increase  in  size  for 
coarser  levels.  As  a  result.  Algorithm  C  is  from  3-7  times  faster  than  Algorithm  D,  for  an 
equivalent  number  of  points  n  and  order  of  convergence  p. 


Chapter  4 

Generalizations  and  Conclusions 


4.1  Generalizations 

In  Chapter  3,  we  decomposed  a  one- dimensional  integral  operator  P  into  four  operators 
PaAi  Pab,  Pbai  Pbb,  and  then  constructed  low  rank  factorizations  of  Pab  and  Pba- 
The  specific  factorization  of  Pab  and  Pba  in  Chapter  3  involved  subdividing  each  of  the 
intervals  A  and  B  into  a  number  of  smaller  subintervals,  and  decomposing  Pab  and  Pba  into 
a  number  of  operators  acting  on  these  smaller  subintervals.  Each  of  the  smaller  operators 
either  acted  on  subintervals  which  were  well-separated  from  each  other  (the  operator  thus 
can  be  approximated  by  a  low-order  Chebyshev  polynomial),  or  acted  on  subintervals  which, 
when  discretized,  contained  few  points  (the  discretization  of  the  operator  thus  being  of  low 
dimension).  Such  factorizations  can  also  be  applied  to  integral  operators  corresponding  to 
integral  equations  on  a  curve,  and  to  operators  corresponding  to  integral  equations  of  two 
and  three  dimensions. 

Two  problems  arise  when  this  method  is  applied  to  integral  equations  on  curves.  First, 
it  is  more  difficult  to  subdivide  A  and  B  into  a  number  of  smaller  well-separated  intervals. 
For  example,  suppose  that  the  curve  is  described  by  a  polynomial,  and  suppose  further  that 
one  needs  to  determine  the  rank  of  the  interaction  between  two  sections  of  the  curve.  It  is 
not  clear  how  to  determine  whether  the  two  sections  are  well-separated.  A  second  problem 
is  that  the  choice  of  boundary  separating  subintervals  A  and  B  can  dramatically  affect  the 
ranks  of  the  operators  Pab  a^iid  Pba-  As  an  example,  consider  an  integral  equation  for  an 
ellipse  in  which  one  of  the  axes  of  the  ellipse  is  much  longer  than  the  other.  If  the  short  axis 
is  chosen  to  be  the  boundary  separating  subintervals  A  and  B,  then  there  are  few  points 
in  the  discretization  which  are  close  to  this  boundary;  Pab  and  Pba  wUl  therefore  be  low 
rank  operators.  On  the  other  hand,  if  the  long  axis  is  chosen  to  be  the  boundary  separating 
A  and  B,  then  nearly  all  the  points  in  the  discretization  will  be  close  to  the  boundary,  and 
the  ranks  of  Pab  and  Pba  will  be  quite  high. 

When  integral  operators  of  two  and  three  dimensions  are  considered,  the  rank  of  the 
operators  Pab  and  Pba  is  largely  determined  by  the  number  of  points  on  the  boundary 
separating  A  and  B.  For  an  iVx  A  discretization  in  two  dimensions,  the  boundary  separating 
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A  and  S  is  a  line  containing  N  points,  and  the  rank  of  Pj^  and  Pba  is  iVlog  Thus  a 
direct  solver  for  a  two  dimensional  integral  equation  would  require  order  0(iV®)  arithmetic 
operations,  and  would  be  of  considerable  practical  interest.  On  the  other  hand,  for  an 
N  X  N  y.  N  discretization  in  three  dimensions,  the  boundary  separating  A  and  5  is  a  square 
containing  N  X  N  points,  and  the  rank  of  Pab  and  Pba  for  this  operator  is  jV^logiV.  A 
direct  solver  for  a  three  dimensional  integral  equation  would  require  order  0(iV®)  arithmetic 
operations,  and  would  be  prohibitively  expensive. 

4.2  Conclusions 

Algorithms  have  been  presented  for  the  solution  of  two-point  boundary  value  problems  for 
ordinary  differential  equations,  and  for  the  solution  of  one-dimensional  first  and  second 
kind  integral  equations  of  potential  theory.  All  algorithms  have  CPU  time  requirements 
proportional  to  N  with  N  the  number  of  nodes  in  the  discretization,  and  p  the  desired 
order  of  convergence.  In  addition,  the  time  requirements  axe  insensitive  to  the  condition 
number  of  the  discretized  linear  system.  The  methods  permit  the  use  of  schemes  with 
extremely  high  orders  of  convergence,  and  are  quite  insensitive  to  end-point  singularities. 
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